COMPLEX  THREE-DIMENSIONAL  TURBULENT  ELOW  PHENOMENA 

OF  HIGH-ALPHA  PROJECTILE  AERODYNAMICS 


GUANG-WU  CHEN 


A DISSERTATION  PRESENTED  TO  THE  GRADUATE  SCHOOL 
OF  THE  UNIVERSITY  OF  FLORIDA  IN  PARTIAL  FULFILLMENT 
OF  THE  REQUIREMENTS  FOR  THE  DEGREE  OF 

DOCTOR  OF  PHILOSOPHY 

UNIVERSITY  OF  FLORIDA 


1995 


ACKNOWLEDGEMENTS 


I wish  to  express  my  sincere  appreciation  to  Dr.  Chen-Chi  Hsu  for  his  guidance 
and  encouragements  throughout  my  research  and  in  working  on  the  dissertation 
manuscript.  I would  also  like  to  extend  my  thanks  to  other  supervisory  committee 
members,  Dr.  W.  Shyy,  Dr.  D.  M.  Mikolaitis,  Dr.  U.  H.  Kurzweg  and  Dr.  J.  Klausner, 
for  their  interest  and  suggestions  in  this  research. 

Special  thanks  are  extended  to  Dr.  Shu-Cheng  Yang,  Dr.  Chau-Lin  Lee,  Dr. 
Shin-Jye  Liang,  Dr.  Edwin  Blosch  and  Mr.  Guobao  Guo  for  their  kind  assistance  since 
I first  joined  the  CED  family  here  in  the  department.  I would  also  like  to  express  my 
gratitude  to  my  family  for  their  love  and  support  over  the  years.  Last  but  not  least  are 
the  people  who  have  helped  me  many  times  and  made  things  possible. 

This  study  was  partially  supported  by  Elorida  State  University  and  by  Pittsburgh 
Supercomputing  Center  through  the  allocation  of  supercomputer  resources. 


11 


TABLE  OF  CONTENTS 


Page 

ACKNOWLEDGEMENTS  ii 

ABSTRACT  vi 

CHAPTERS 

I INTRODUCTION  1 

High-a  Aerodynamics 3 

Role  of  CFD 6 

High-a  Projectile  Aerodynamics 7 

Crossflow  Separation  7 

Three-Dimensional  Flow  Separation  8 

Topology  of  Three-Dimensional  Flow  Separation  9 

High-a  Flow  Regimes 11 

Literature  Survey 11 

Research  Objectives  19 

II  GOVERNING  EQUATIONS  21 

Navier-Stokes  Equations  and  Equations  of  State 21 

Continuity  Equation 22 

Momentum  Equations 22 

Energy  Equation 23 

Governing  Equations  in  Vector  Form 24 

Equations  of  State 25 

Non-dimensional  Navier-Stokes  Equations 26 

Transformed  Navier-Stokes  Equations 28 

Thin-Layer  Navier-Stokes  Equations 35 

Turbulence  Model 36 

Modified  Baldwin-Lomax  Turbulence  Model  41 

III  FLOW  PROBLEM 44 


111 


Flow  Conditions  and  Geometrical  Model 44 

Domain  Discretization 45 

Planar  Grid  Generation 46 

Three-Dimensional  Grid  Generation  46 

Boundary  Conditions 47 

Solid  Surface  Boundary  Condition  47 

Outer  Domain  Boundary  Condition 49 

Circumferential  Boundary  Condition 49 

Singular  Plane  J = 1 49 

Downstream  Boundary  Condition 50 

Slow  Start 51 

IV  THIN-LAYER  NAVIER-STOKES  COMPUTATION 55 

Thin-Layer  Navier-Stokes  TVD  Code 56 

Solution  Method  and  Application  59 

Algorithm  and  Flow  Chart 60 

Residual 62 

Pressure  and  Shear  Stress  Coefficients 62 

Aerodynamic  Force  66 

Mach  Number  67 

Vorticity 68 

Streamtraces  and  Skin-Friction  Lines 68 

Evaluation  of  the  Jacobian  J 70 

Integration  Time  Step  At  70 

V EFFECTS  OF  TURBULENCE  MODELS  73 

Impacts  of  Excessive  Eddy  Viscosity  74 

On  the  Crossflow  Plane  J=54 74 

On  the  Surface  Flow  Properties 78 

On  the  Development  of  Crossflow  Separation 83 

Implication  of  the  Cut-off  Constant 87 

VI  TURBULENT  VORTICAL  FLOW  PHENOMENA 114 

Onset  of  Cross-sectional  Vortical  Flow  Structure 115 

Onset  and  Development  of  Surface  Flow  Patterns 118 

Three-Dimensional  Streamtraces 122 

Vortex  Shedding  Phenomena 123 

VII  SUMMARY  AND  CONCLUDING  REMARKS 163 

APPENDIX:  JACOBIAN  MATRICES  168 


IV 


REFERENCES 170 

BIOGRAPHICAL  SKETCH 175 


V 


Abstract  of  Dissertation  Presented  to  the  Graduate  School 
of  the  University  of  Florida  in  Partial  Fulfillment  of  the 
Requirements  for  the  Degree  of  Doctor  of  Philosophy 

COMPLEX  THREE-DIMENSIONAL  TURBULENT  FLOW  PHENOMENA 

OF  HIGH-ALPHA  PROJECTILE  AERODYNAMICS 

By 

Guang-Wu  Chen 
August  1995 

Chairperson:  Dr.  Chen-Chi  Hsu 

Major  Department:  Aerospace  Engineering,  Mechanics  and  Engineering  Science 

An  impulsively  started  turbulent  flow  of  Mach  0.96  past  a Secant-Ogive-Cylinder- 
Boattail  projectile  model  with  sting  at  60-degree  angle  of  attack  is  considered  for 
investigating  the  complex  three-dimensional  vortical  flow  physics  and  phenomena.  The 
governing  equations  used  for  the  vortical  flow  simulation  are  the  Reynolds-averaged  thin- 
layer  Navier-Stokes  equations  implemented  with  a modified  Baldwin-Lomax  algebraic 
eddy  viscosity  model  for  ideal  gases.  The  turbulence  model  has  been  shown  to  provide 
an  acceptably  accurate  surface  pressure  distribution  in  the  crossflow  separation  region  by 
Degani  and  Schiff.  The  implication  of  the  turbulence  model  for  a proper  description  of 
flow  solutions  away  from  the  boundary-layer  region  has  been  intensively  investigated  in 
this  study.  The  comparison  of  the  simulated  flow  characteristics  against  those  obtained 
from  employing  the  original  Baldwin-Lomax  model  clearly  shows  that  the  modified 
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turbulence  model  employed  in  this  study  can  give  at  least  qualitatively  correct  vortical 
flow  physics  and  phenomena. 

In  this  study  the  numerical  simulation  has  revealed  a number  of  interesting  and 
important  fundamental  characteristics  of  complex  vortical  flow.  The  onset  times  for 
primary  and  secondary  crossflow  separations  are  found  to  be  linearly  dependent  on  the 
diameter  of  the  cross-section.  The  evolution  of  the  crossflow  vortex  patterns  on  the 
cylinder,  boattail  and  sting  sections  of  the  projectile  shows  a transition  from  a stable  node 
to  an  unstable  spiral  through  the  appearance  of  a stable  limiting  cycle.  The  unstable 
spiral  vortex  on  the  ogive;  however,  is  formed  through  an  entirely  different  mechanism, 
consisting  of  a two-saddle-point  flow  structure.  It  is  also  shown  that  the  development  of 
primary  vortices  has  induced  a high  speed  axial  jet  stream  in  the  leeward  side  of  the 
ogive  and  cylinder  sections.  Moreover,  the  alternating  vortex  shedding  is  simulated  on 
the  projectile  sting  and  the  shedding  frequency  in  terms  of  the  Strouhal  number  ranges 
from  0.179  to  0.191  on  different  cross-sections. 
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CHAPTER  I 
INTRODUCTION 


In  the  past  few  decades,  aircraft  designers  have  been  continuously  creating  flight 
vehicles  of  improved  performance,  and  with  broader  flight  envelope.  Modern  jet  fighters 
of  the  1990s  can  fly  much  faster  and  higher  than  their  counterparts  in  World  War  II. 
The  advancement  of  avionics  and  weapon  systems  allows  them  to  carry  out  diverse 
missions.  Stealth  technology  can  even  make  planes  invisible  to  enemy  radar.  No  matter 
how  far  and  how  high  a modern  Jet  fighter  can  go,  the  most  vital  part  comes  back  to  a 
simple  matter:  the  traditional  "dog  fight"  (close-in  combat)  in  the  range  of  Mach  number 

0.8  to  1.2.  The  modern  fighter  must  be  highly  maneuverable  to  counteract  its  foe  in  the 
sky. 

In  1989  a Russian  demonstrator  Su-27  performed  a famous  "Cobra  Maneuver" 
in  an  air  show  and  stunned  the  observers  [1].  The  airplane’s  nose  pitched  up  beyond  the 
vertical  plane  but  the  airplane  was  able  to  maintain  a level  path  and  recover  its  normal 
flight.  The  airspeed  dropped  almost  200  knots  («230  mph)  within  a few  seconds.  Such 
strong  agility  can  give  the  plane  a great  advantage  in  a dog  fight,  thanks  to  its  excellent 
high-angle-of-attack  performance.  How  the  Su-27  can  maintain  control  at  such  an 

extremely  high  angle  of  attack  has  created  much  interest  for  western  designers  in  recent 
years . 
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All  aerodynamic  engineers  are  familiar  with  one  simple  rule:  for  a given  flight 
condition,  the  airfoil  lift  increases  when  the  angle  of  attack  a goes  up.  But  there  exists 
a critical  a.  where  the  lift  reaches  its  maximum.  Above  that  critical  a,  the  flow  on  top 
of  the  airfoil  separates  and  causes  a pressure  back  up.  With  the  absence  of  a low 
pressure  suction  on  the  upper  airfoil,  the  lift  drops  suddenly  and  leads  to  a stall.  For  an 
airplane  the  lift  can  also  be  contributed  from  other  components  such  as  the  nose  and 
fuselage.  When  a is  large,  the  incoming  flow  separates  from  both  sides  of  the  slender 
nose/fuselage  configuration  and  forms  a pair  of  highly  counter-rotating  vortices  on  the 
leeward  side.  The  low  pressure  region  developed  by  the  vortical  flow  can  provide 
additional  lift  which  improves  the  airplane’s  maneuverability  and  agility.  Vortical  flow 
does  have  some  adverse  effects  on  the  flight  vehicle,  such  as  lateral  instability  due  to  the 
asymmetric  vortex  pattern.  A fundamental  knowledge  of  the  vortex  dynamics  is 
therefore  desirable  to  fully  take  advantage  of  the  vortical  flow,  while  in  the  meantime, 
minimizing  the  possible  adverse  effects.  In  current  study,  the  general  high-a  vortical 
flow  phenomena  are  investigated  through  a numerical  simulation,  in  particular,  a 
transonic  turbulent  flow  past  a projectile  configuration.  In  the  following,  some  general 
aspects  of  the  high-a  vortical  flow  problem  toward  an  airplane  are  first  addressed.  A 
detailed  background  of  high-a  projectile  aerodynamics,  including  three-dimensional  flow 
separation,  flow  topology,  and  a literature  survey  is  presented,  which  is  followed  by  the 
goal  of  this  study. 
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High-o;  Aerodynamics 

Aircraft  lateral  (or  directional)  stability  and  control  characteristics  at  high  angle 
of  attack  can  be  significantly  altered  by  the  presence  of  vortices  developed  from  the 
aircraft  nose  and  forward  fuselage.  Because  forebody  vortices  at  high  a’s  are  often 
asymmetric  relative  to  the  angle-of-attack  plane,  their  action  on  body,  wing,  and  tail 
surfaces  can  result  in  large  yawing  and  rolling  moments.  Such  effects  may  cause 
conventional  stability  and  control  surfaces  (rudders,  ailerons  etc.)  to  be  less  effective  at 
high  angle  of  attack.  Another  example  is  the  tail  buffet  phenomenon  resulting  from  the 
tail’s  interaction  with  the  shedded  forebody  vortices,  exciting  the  natural  frequencies  of 
the  fin(s)  structure.  The  interactions  between  asymmetric  forebody  vortices  with  other 
airplane  surfaces  at  sideslip  conditions  can  also  lead  to  the  wing  rock  phenomenon— a 
sustained  periodic  oscillation  in  pitch  and  roll.  These  oscillations  may  reduce  the 
structural  service  life,  even  causing  a sudden  structural  failure  during  flight.  The  vortical 
flow  structure  can  also  breakdown  such  that  the  lift  force  generated  by  the  vortices  is 
suddenly  lost,  causing  a possible  aircraft  stall. 

Control  of  forebody  vortex  asymmetry  is  therefore  of  interest  to  minimize  the 
adverse  effects  and  to  use  the  associated  yawing  and  rolling  moments  to  enhance  actively 
the  aircraft’s  high-angle-of-attack  performance.  Nose  strakes,  chines,  or  leading-edge- 
extension  (LEX)  are  found  on  several  aircrafts.  These  devices  provide  additional  lift  due 
to  the  vortical  flow  they  develop  at  moderate-to-high  angles  of  attack.  These  are  useful 
for  minimizing  forebody  vortex  asymmetry;  however,  the  prediction  and  control  of  the 
vortical  flow  and  mutual  interactions  of  the  vortices  are  not  well  known.  Asymmetric 
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nose  strakes  and  asymmetric  air  jet  injection  can  induce  a particular  orientation  and 
degree  of  forebody  vortex  asymmetry,  but  such  behavior  and  the  resultant  effects  on  the 
complete  aircraft  again  are  not  yet  understood  completely.  The  combined  effect  of  the 
LEX  or  wing-body  strake  vortices,  as  well  as  the  forebody  vortices  on  the  vehicle 
aerodynamics,  must  be  considered  to  avoid  the  flight  stability  control  and  other  problems 
such  as  tail  buffet  or  wing  rock.  The  appearance  of  shock  waves  in  the  high  speed  range 
results  in  some  phenomena  associated  with  shock/boundary-layer,  shock/vortex 
interactions  that  can  also  complicate  the  flowfield. 

Recently,  there  has  been  growing  interest  from  government  agencies  and  industry 
to  tackle  those  issues.  At  NASA  Ames  Research  Center,  extensive  experimental  and 
computational  investigations  have  been  conducted  on  the  High  Alpha  Research  Vehicle 
(HARV)  to  determine  the  effects  of  three-dimensional  separated  vortical  flow  on  the 
aircraft  and  means  of  utilizing  these  effects  to  the  best  possible  advantages  [2].  The  F/A- 
18  was  selected  as  the  HARV  base  configuration  because  of  its  exceptional  high-angle-of- 
attack  capability.  The  combination  of  its  complex  geometry  and  the  flow  physics 
associated  with  the  vortical  flow  at  high  a makes  both  the  experimental  and 
computational  processes  a great  challenge  to  the  engineering  society.  Significant  results 
have  been  obtained  and  reported  in  several  conference  papers,  for  example.  Refs.  3-5. 

In  addition  to  the  F-18  HARV  at  NASA,  the  forward-swept-wing  X-29 
experimental  aircraft  is  also  currently  under  consideration  for  flight  testing  of  high  a 
forebody  vortex  control  concepts  at  Dryden  Flight  Research  Center  [6,7].  This 
technology  demonstrator  was  designed  to  be  maneuverable  at  an  angle  of  attack  up  to 
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70°.  The  issues  of  the  study  include  (1)  the  ability  to  predict  asymmetric  vortical 
forebody  flow,  (2)  the  use  of  asymmetric  forebody  jets  to  induce  formation  of 
asymmetric  vortices,  and  (3)  the  use  of  asymmetric  forebody  Jets  to  reduce  vortex 
asymmetry.  The  X-31  Enhanced  Fighter  Maneuverability  (EFM)  demonstrator  is  another 
successful  example  of  high-a  aerodynamics  [8].  To  maintain  effective  control  during  a 
post-stall  flight  at  extremely  high  angles  of  attack,  the  "thrust  vector  vanes"  encircling 
the  engine  exhaust  nozzle  and  a pair  of  "floating"  canards  were  used  on  a delta-wing 
fighter  configuration.  Together  with  a sophisticated  digital  flight  control  system,  the  X- 
31  was  found  to  be  able  to  sustain  a controlled  flight  at  a =70°  for  53  seconds 
(September  18,  1992).  Later  it  achieved  the  first  "Herbst"  maneuver— a 180°  reversal 
of  the  flight  path  carried  out  in  a deep  stall  condition  at  the  speed  of  90  knots  (*  104 
mph)  and  a minimum  50°  angle  of  attack,  after  deceleration  from  a high-speed  entry, 
with  subsequent  high  speed  exit  in  the  opposite  direction  (April  29,  1993).  Such  super- 
agility can  surely  give  a fighter  aircraft  great  advantages  in  close  air  combat. 

Erickson  [9]  discussed  both  civil  and  military  airplane  high-a  aerodynamics  with 
a large  variety  of  configurations.  In  a much  broader  sense  he  defined  the  high-a  as  any 
situation  where  an  airplane  or  any  of  its  components  encounters  an  off-design,  high-lift 
flight,  including  attached  flows,  uncontrolled  separated  flows,  and  organized  flow 
separation  in  the  form  of  vortices,  shock  waves,  or  combinations  of  these  flow 
conditions.  For  example,  the  high-a  flow  phenomena  can  be  seen  in  the  propeller  tip 
vortex  system  on  a light  general  aviation  airplane,  blade  tip  vortex  system  on  a 
rotorcraft,  and  wing-mounted  engine  nacelle  and  pylon  assembly  of  a transport  airplane. 
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A number  of  flow  control  concepts  such  as  airframe  shaping  or  actively  controlled 
mechanical/pneumatic  devices  which  effectively  improve  the  quality  of  high-a 
aerodynamics  were  reviewed. 


Role  of  CFD 

For  decades,  wind-tunnel  or  flight  tests  have  been  major  sources  of  engineering- 
level  data  to  aircraft  designers.  The  types  of  data  available  from  wind-tunnel  tests  are 
usually  limited  by  the  cost  and  instrumentation  restrictions.  The  measured  data  often 
lack  resolution  and  might  be  affected  by  the  existence  of  the  measuring  devices  sitting 
in  the  flowfield.  The  flight  test  costs  too  much  and  gives  even  less  data  because  of  larger 
instrumentation  restrictions.  For  some  most  interesting  cases  such  as  the  high-a 
maneuver,  it  is  simply  too  dangerous  to  fly.  In  addition  to  the  field  or  flight  test, 
computational  fluid  dynamics  (CFD)  is  another  approach  which  can  provide  valuable 
information  in  the  design  process.  A qualitative  description  to  capture  the  dominant  flow 
phenomena  is  often  more  important  than  an  absolutely  accurate  prediction  in  the 
preliminary  stage.  Once  we  are  on  a correct  track,  the  code  can  be  improved  through 
the  numerical  validation  which  requires  comparison  with  carefully  conducted  experiments 
[10].  The  validation  of  a CFD  code  with  experiments  can  lead  to  a general  confidence 
in  a given  methodology  to  capture  the  flow  phenomena  of  interest.  On  the  other  hand, 
the  more  detailed  information  obtained  by  a flow  simulation  can  help  to  interpret  and 
clarify  some  experimental  results.  With  the  advent  of  supercomputers  and  the 
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advancement  of  numerical  techniques  for  nonlinear  problems,  CFD  now  can  provide 

more  comprehensive  flow  physics  and  phenomena  than  the  conventional  experiment  does 

if  the  mathematical  model  solved  is  a suitable  one  for  the  complex  flow  problem. 

The  creative  use  of  CFD  methods  can  address  the  flow  physics  of  design 
problems  and  configuration  concepts  for  which  no  design  experiences,  data  base, 
or  test  techniques  exist.  This  approach  can  both  yield  critical  design  information 
and  identify  areas  of  high  potential  for  future  methodology  development  and 
applications  programs.  [6,  pp.  818] 


High-g  Projectile  Aerodynamics 

The  study  of  a viscous  flow  past  a projectile  geometry  (a  slender  body  of 
revolution)  at  high  angles  of  attack  is  essential  to  the  understanding  of  the  high-g  aircraft 
aerodynamics.  The  complex  vortical  flow  phenomena,  such  as  the  three-dimensional  flow 
separation  and  the  development  of  the  symmetric/asymmetric  vortical  flow  in  the  wake 
region  which  arise  on  a highly  maneuverable,  geometrically  complicated  flight  vehicle 
at  high  g,  can  also  be  observed  regardless  the  simplicity  of  the  projectile  shape.  Besides, 
it  is  much  easier  to  apply  CFD  techniques  such  as  grid  generation  and  graphical 
representation  to  a simple  geometry  like  this.  The  general  characteristics  of  a vortical 
flow  behind  an  inclined  projectile  body  are  described  as  follows. 

Crossflow  Separation 

A schematic  vortical  flow  structure  above  a projectile  body  is  shown  in  Figure 
1.1  for  a crossflow  plane.  The  outer  flow  approaching  the  windward  plane  of  symmetry 
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turns  and  flows  outward  along  the  body  from  the  windward  side  to  the  leeward  side.  A 
well-attached  boundary  layer  is  developed  until  it  separates  from  the  body  at  a primary 
separation  point  Sj;  accordingly,  a pair  of  counter-rotating  primary  vortices  on  the 
leeward  side  of  the  body  is  formed.  The  primary  vortices  induce  a flow  downward  to 
the  body  surface  which  attaches  onto  the  leeward  plane  of  symmetry  at  A,.  The  induced 
downflow  turns  outward  along  the  body  surface  from  the  leeward  attachment  point  to  the 
primary  separation  Sj.  The  adverse  circumferential  pressure  gradient  under  the  primary 
vortices,  if  strong  enough,  can  cause  the  induced  flow  to  separate  at  a secondary 
separation  point  S2.  The  local  fluid  rolls  up  behind  the  body  to  form  the  secondary 
vortex  structure.  Another  point  of  attachment  A2  is  located  on  the  body  surface  between 
the  primary  and  secondary  separation  points.  Figure  1.1  shows  that  the  separating  flow 
forms  an  unstable  vortex  structure,  i.e.,  the  fluid  goes  away  from  the  vortex  center.  It 
is  mentioned,  however,  that  the  leeward  side  crossflow  structure  can  consist  of  either 
stable/unstable  vortices  or  limiting  cycles  depending  on  the  various  three-dimensional 
flow  development  [11]. 

Three-Dimensional  Flow  Separation 

In  the  two-dimensional  flow  case,  for  example  the  flow  past  a cylinder,  there 
exists  only  one  shear  stress  component  along  the  circumferential  direction.  The 
separation  or  reattachment  then  can  be  easily  located  as  the  surface  points  at  which  the 
shear  stress  becomes  zero.  For  three-dimensional  flow,  in  addition  to  the  crossflow 
shear  stress  component,  an  axial  component  can  carry  the  flow  downstream  to  form  a 
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line  of  separation  which  is  an  asymptotic  convergence  of  adjacent  surface  skin-friction 
lines.  The  combined  action  of  the  circumferential  and  axial  components  on  the  body 
surface  may  be  seen  in  Figure  1.2,  which  shows  the  simulated  surface  skin-friction  lines 
on  part  of  an  unwrapped  projectile  body.  The  axial  component  in  this  graph  is  to  the 
right  hand  side.  The  convergence  of  surface  skin-friction  lines  toward  the  primary  and 
secondary  separation  lines,  Sj  and  S2,  and  the  circumferential  divergence  of  the  skin- 
friction  lines  from  the  attachment  lines.  A,  and  A2,  are  clearly  visible.  It  should  be 
clarified  that  the  flow  does  not  have  to  "leave"  the  surface  into  the  wake  region  on  the 
separation  lines  with  the  circumferential  shear  stress  component  zero.  In  two- 
dimensional  flow  the  change  sign  of  shear  stress  simply  means  the  fluid  must  be  elevated 
and  form  a recirculation  bubble  above  the  surface.  In  the  three-dimensional  flow  case, 
the  fluid  is  allowed  to  escape  in  the  downstream  direction  to  form  the  skin-friction 
convergence  lines  as  shown  here.  This  makes  the  study  of  the  surface  flow  topology 
important  in  viewing  and  understanding  a three-dimensional  flow  separation. 

Topology  of  Three-Dimensional  Flow  Separation 

One  of  the  approaches  to  study  three-dimensional  flow  separation  is  to  make  use 
of  the  topology  with  the  aid  of  critical-point  theory  [12].  Also  known  as  the  "phase- 
plane"  theory,  the  critical-point  theory  was  first  used  in  solving  autonomous  ordinary 
differential  equations.  The  detailed  classifications  of  various  types  of  separated  flow, 
nomenclature,  and  topological  terminology  were  described  in  a review  paper  by  Chapman 
and  Yates  [13].  In  the  following,  only  materials  directly  related  to  the  current  vortical 
flow  problem  are  presented. 
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Mathematically  critical  points  are  often  referred  to  as  singular  points.  A singular 
point  in  a velocity  field  implies  that  the  magnitude  of  the  velocity  vector  at  the  point  is 
zero,  i.e.,  the  slope  of  particle  trajectories  is  indeterminate.  For  a viscous  flow  on  the 
solid  surface  where  the  velocity  vanishes  due  to  no-slip,  the  two-component  shear  stress 
vector  is  used  instead  of  a 3-D  velocity  vector.  According  to  their  behavior,  critical 
points  are  known  as  either  nodal  points  or  saddle  points.  A nodal  point  is  characterized 
by  an  infinite  number  of  trajectories  entering  (stable  node)  or  leaving  (unstable  node)  the 
point.  It  can  be  further  classified  as  either  a regular  node  or  a spiral  node.  A saddle 
point  has  only  two  lines  called  separatrixes  intersect  at  the  point,  one  is  going  toward  the 
point,  and  the  other  one  away  from  it.  There  are  rules  that  can  be  followed  on  3-D  body 
surface  or  crossflow  planes  regarding  the  number  of  critical  points  [14].  On  the  surface 
of  a closed  3-D  body,  the  summation  rule  states  that 

'£N~Y^S  = 2 (1-1) 

where  N denotes  the  isolated  node  and  S the  saddle  point.  On  a crossflow  plane,  the 
summation  rule  becomes 

N + -YN'  - 

where  N’  and  S’  denote  half-node  and  half-saddle  that  occur  on  a boundary.  There  is 
also  a local  summation  rule  for  part  of  surface  skin-friction  line  portrait  which  can  be 
isolated  and  enclosed  within  a smooth  curve.  This  rule  says  that,  within  the  isolated 
region,  a saddle  point  must  be  added  or  deleted  when  each  node  is  added  or  deleted. 
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The  interpretation  of  surface  flow  topology  can  be  applied  to  the  3-D  flowfield. 
For  an  unstable  node  on  the  solid  surface,  the  trajectories  are  all  pointing  away  from  the 
node  which  suggests  the  existence  of  a source  of  fluid  flow  on  the  wall.  This  implies  that 
the  flow  attaches  to  the  node  and  spreads  out  along  the  surface.  Such  a node  is  called 
a node  of  attachment  N^.  A node  of  separation  Ns,  therefore,  accounts  for  a stable  node 
which  is  a sink-like  point  on  the  wall.  Similarly  the  saddle  of  separation  Sg  and  saddle 
of  attachment  are  defined  as  detailed  in  Ref.  13. 

High-g  Flow  Regimes 

The  leeward  side  flow  pattern  on  an  inclined  ogive-cylinder  projectile  can  be 
usually  categorized  into  three  regimes.  For  a is  small,  say  less  than  30°,  a symmetric 
vortical  flow  pattern  is  observed  with  the  crossflow  topology  being  similar  to  Figure  1.1. 
In  general  the  vortex  size  grows  with  the  longitudinal  distance  from  the  nose.  For  a is 
30°  —60°,  the  steady  asymmetric  vortex  pattern  is  observed.  For  g is  60°  and  above, 
vortex  shedding  becomes  essential,  with  an  oblique  shedding  near  the  nose  and  parallel 
shedding  far  downstream.  With  g approaching  90°,  the  parallel  shedding  turns  out  to 
be  dominant  over  the  oblique  shedding  except  very  close  to  the  nose  tip. 


Literature  Survey 

Over  the  past  decade  there  were  many  research  efforts  into  high-g  projectile 
aerodynamics.  These  works  included  experimental  or  computational  approaches  on 
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laminar  or  turbulent  flows  ranging  from  subsonic  to  supersonic  regime.  A number  of 
selected  works  are  discussed  in  the  following. 

In  1986,  Degani  and  Schiff  [15]  used  a thin-layer  parabolic  Navier-Stokes  method 
to  compute  turbulent  supersonic  flows  around  various  pointed  bodies  at  large  incidence 
up  to  22.75°.  Due  to  the  significant  crossflow  separation,  they  first  tested  whether 
neglecting  the  circumferential  viscous  terms  in  the  thin-layer  approximation  could 
degrade  the  solution  accuracy.  Adding  back  those  terms  caused  an  1 1 % increase  in  CPU 
time,  but  produced  less  than  0.5  % change  in  the  computed  solution.  Practically,  the 
circumferential  grid  resolutions  are  not  sufficient  to  resolve  the  circumferential  viscous 
terms.  Secondly  they  proposed  a modification  of  the  well-known  Baldwin-Lomax 
turbulence  model  [16]  due  to  the  large  crossflow  separation,  which  effectively  reduced 
the  over-predicted  eddy  viscosity.  The  computed  surface  pressure  distribution  matched 
experimental  data  much  better  than  the  one  predicted  by  the  original  model.  Also  in 
some  computed  cases,  the  secondary/tertiary  vortices  were  clearly  visible  which  were 
absent  from  the  original  model  computation.  They  also  performed  the  computation  by 
different  circumferential  discretization  of  5°  and  2.5°,  respectively.  The  results  showed 
that  a grid  size  of  2.5°  is  required  to  resolve  the  small-scale  secondary  vortex  structure. 
The  details  of  the  Baldwin-Lomax  model  and  modification  are  to  be  discussed  in  the  next 
chapter. 

A further  investigation  on  the  effect  of  the  modified  Baldwin-Lomax  turbulence 
model  and  circumferential  grid  resolution  to  the  vortical  flow  simulation  was  conducted 
by  Schiff,  Degani  and  Cummings  [17].  In  the  study  they  utilized  a time-marching  thin- 
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layer  Navier-Stokes  code  to  compute  Mach  number  0.2  flow  past  a projectile  with  3.5D 
tangent-ogive  forebody  and  7.0D  cylindrical  afterbody  at  a =20°,  25°  and  30°.  The 
flow  was  assumed  to  be  symmetric  and  only  half  of  the  domain  was  solved  with  the 
circumferential  discretizations  of  A0=3°  and  1.5°.  The  computed  surface  flow  patterns 
showed  the  primary  crossflow  separation  line  predicted  with  the  modified  model  was 
located  more  upstream  than  the  one  predicted  by  the  original  model.  Furthermore,  the 
modified  model  computation  showed  the  development  of  a tertiary  separation  line  near 
the  rear  body.  On  observing  the  crossflow  vortical  flow  pattern,  the  modified  model 
computation  showed  a larger,  more  complex  structure.  By  using  a grid  resolution  of 
A9—3°  at  o!=20°,  it  was  found  that  the  code  had  trouble  in  differentiating  between  the 
primary  and  secondary  separation  lines  on  the  ogive  part.  Reducing  the  grid  spacing  by 
half  cured  the  problem.  However,  the  requirement  of  increasing  circumferential  grid 
resolution  was  significantly  reduced  as  the  flow  physics  were  altered  at  larger  angles  of 
attack.  At  a =25°  and  30°,  the  predicted  primary  and  secondary  separation  lines  were 
further  apart  due  to  a stronger  circumferential  adverse  pressure  gradient.  Therefore  a 
grid  spacing  of  A6=3°  was  sufficient  to  resolve  the  crossflow  structure.  Their 
experience  suggested  that  at  least  4 or  5 circumferential  points  are  needed  between 
primary  and  secondary  separation  lines  in  order  to  realistically  model  the  flowfield. 

As  has  been  mentioned  previously,  the  leeward  side  flow  becomes  very  unstable 
at  large  angles  of  attack.  Degani  and  Zilliac  [18]  performed  the  surface  pressure  and 
wake  velocity  measurements  in  a low-turbulence  wind  tunnel  past  an  ogive-cylinder  body. 
The  flow  was  laminar,  with  a Reynolds  number  of  26,000.  The  flow  was  found  to  be 
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stable  with  a less  than  15°  to  20°.  The  flowfield  at  a=30°  was  found  to  be  highly 
unsteady  by  observing  some  high  frequency  signals  (5000  — 6000  Hz)  with  low 
amplitude,  which  was  associated  with  the  small-scale,  three-dimensional  vortices  in  the 
shear  layer.  As  a increased  from  30°  to  60°,  a low  frequency  component  («200  Hz) 
started  to  develop  and  became  identifiable  at  a =60°  — 85°  near  the  rear  of  the  cylindrical 
afterbody.  The  measured  frequency  was  130— 176  Hz  which  was  found  to  be  related  to 
the  von  Karman  vortex  shedding.  The  non-dimensional  frequency  or  Strouhal  number 
was  found  to  be  0.21sin(a)  where  0.21  is  the  classical  Strouhal  number  of  the  von 
Karman  shedding  measured  on  a two-dimensional  cylinder  [19]. 

At  a moderate  to  high  a before  the  vortex-shedding  exists,  the  small  unsteadiness 
is  able  to  trigger  the  asymmetry  and  maintains  a "steady"  asymmetric  vortex  pattern. 
How  a symmetric  body  can  produce  a non-symmetric  flow  pattern  has  caused  great 
interest.  Degani  and  Schiff  [19]  performed  a thin-layer  Navier-Stokes  laminar  vortical 
flow  simulation  of  Mach  0.2  flow  past  an  ogive-cylinder  body  at  40°  angle  of  attack. 
The  circumferential  discretization  had  A6=3°  for  the  whole  domain.  The  computed  flow 
remained  symmetric,  until  a symmetry-breaking  perturbation  was  applied  near  the  nose 
to  induce  the  asymmetric  pattern.  When  the  perturbation  was  removed,  the  asymmetric 
solution  returned  to  symmetric.  A space-fixed,  time-invariant  perturbation  was  needed 
to  reach  and  maintain  an  asymmetric  flow  pattern.  A stronger  perturbation  was  found 
to  result  in  a greater  asymmetry,  and  so  does  a perturbation  inserted  at  an  upstream 
location  closer  to  the  nose  tip.  Two  pairs  of  vortices  were  found  which  were  consistent 
with  what  experimentally  had  been  observed  [21]. 
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Degani  [22]  continued  the  computation  for  a =60°  and  80°  by  inserting  a small 
transient  perturbation  near  the  tip.  Unlike  the  a=40°  case,  the  triggered  asymmetric 
flow  pattern  did  not  return  to  symmetry  when  the  perturbation  was  removed.  The  flow 
remained  highly  asymmetric  instantaneously  with  the  computed  side  force  history 
showing  oscillations  around  a symmetric  mean.  Once  the  perturbation  stayed,  the 
fluctuating  flow  resulted  in  an  asymmetric  mean  side  force.  The  vortiees  were  shed  out 
like  von  Karman  vortex  street.  He  concluded  that  two  kinds  of  instability  can  exist  in 
the  computation.  For  2-D  cylinder  flow  or  3-D  finite-length  body  at  very  large  a (e.g., 
80°),  the  developing  of  von  Karman  vortex  shedding  is  a result  of  absolute  instability. 
A short-time  symmetry-breaking  perturbation  is  required  to  trigger  the  unsteady  vortex 
shedding  phenomenon.  For  the  moderate  to  high  a regime  (a =40°  and  60°),  on  the 
other  hand,  a permanent  disturbance  is  needed  in  order  to  reach  a steady,  asymmetric 
vortical  flow  pattern  as  has  been  observed  in  flow  visualization.  This  is  ealled  a 
convective  instability  of  the  symmetric  flow.  These  two  instabilities  may  coexist  in  a 
certain  range  of  a,  say  60°  to  80°. 

In  fact  geometrical  perturbation  does  exist  in  the  experimental  model.  The 
imperfection  due  to  the  surface  roughness  simply  cannot  be  completely  eliminated. 
Zilliae,  Degani  and  Tobak  [23]  conducted  a series  of  experiments  of  an  ogive-cylinder 
projeetile  to  investigate  the  sensitivity  of  the  leeward  side  flow  structure  to  geometrical 
alteration.  The  ogive  tip  was  made  rotatable  by  a microstepping  motor  and  the  side 
forces  were  recorded  at  different  rotating  angles  0°  -360°,  at  various  angles  of  attack. 
The  results  showed  that  the  variation  of  side  force  with  roll  angle  of  the  nose  was  first 
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found  when  a >30°.  When  a increased,  the  magnitude  of  side  force  variation  with  roll 
angle  grew  from  a continuous  distribution  to  one  gradually  approaching  a square-wave 
distribution.  A nearly  pure  bistable  condition  was  found  at  a =53°,  in  which  the  leeside 
flow  structure  had  only  two  orientations.  At  a =65°  or  larger,  a periodic  vortex 
shedding  was  first  found,  and  the  side  force  variation  with  the  roll  angle  was  becoming 
virtually  zero  like  in  the  low-angles-of-attack  case.  The  study  showed  that  the 
microasymmetry  near  the  nose  acts  as  a permanent  geometrical  perturbation  to  trigger 
the  instability  at  intermediate  angles  of  attack.  On  studying  the  surface  flow  pattern, 
Degani,  Tobak  and  Zilliac  [24]  suggested  that  the  bistable  nature  observed  between 
a =50°  and  65°  is  related  to  the  foci  appearing  near  the  ogive-cylinder  juncture.  In  the 
experiment,  oil  was  attracted  to  the  foci  and  created  thick  bubbles  only  when  a was  up 
to  50°.  As  a further  increased,  the  foci  moved  forward  to  the  tip  and  disappeared  as  a 
approached  70°.  This  phenomenon  was  observed  over  the  testing  Reynolds  number  of 
17,000  — 35,000  in  a laminar  flow  regime. 

The  computational  investigations  mentioned  above  all  assert  that  a symmetry- 
breaking  perturbation  must  be  introduced  in  order  to  achieve  an  asymmetric  solution  as 
had  been  observed  in  experiments  at  moderate  angles  of  attack.  Siclari  and  Marconi  [25] 
performed  the  laminar  Navier-Stokes  simulation  over  a 5°  circular  cone  (<^=5°)  with  the 
assumption  of  conical  flow  solution.  The  testing  case  was  Mach  number  1.8  and 
Reynolds  number  10^  The  grid  used  had  a circumferential  discretization  of  A0=2.25°. 
They  argued  that  the  assumption  of  the  conical  flow  will  not  affect  the  essential  features 
of  the  asymmetric  vortex  flow.  For  a/0 =2,  the  solution  was  able  to  converge  to  a 
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symmetric  flowfield  with  the  residual  dropping  to  machine  zero  and  staying  at  that  level. 
For  a/4>=3  and  up,  the  residual  dropped  almost  10  orders  with  a symmetric  solution, 
then  bounced  back  to  its  original  level  and  then  declined  again.  Finally  the  residual 
reached  machine  zero  and  stayed,  with  the  solution  being  highly  asymmetric.  When  they 
started  the  computation  with  a slightly  asymmetric  initial  condition  (a  2°  sideslip),  it 
converged  monotonically  to  machine  zero,  with  the  identical  asymmetric  solution 
obtained  from  the  symmetric  initial  condition.  They  also  computed  the  resulting  side 
force  and  found  excellent  agreement  with  the  experimental  data.  The  computed  flow 
returned  to  symmetry  when  a/<l)  > 5. 

Siclari  [26]  further  investigated  this  problem  with  various  cone  angles,  Mach 
numbers  and  cross-sectional  shapes.  He  found  that  the  existence  of  the  asymmetric 
solutions  was  affected  by  both  Mach  number  and  cone  angle.  When  the  cone  angle 
increases,  the  asymmetry  declines  very  rapidly.  When  Mach  number  increases,  the 
solution  also  has  a tendency  to  be  symmetric.  He  concluded  that  the  appearance  of  the 
supersonic  crossflow  region  with  the  increasing  free  stream  Mach  number  can  diminish 
the  communication  between  the  opposing  primary  vortical  flow  structures;  therefore,  the 
flowfield  returns  to  symmetry.  Regarding  the  effect  of  the  cross-sectional  shape,  a flatter 
one  (larger  aspect  ratio)  was  found  to  cause  less  asymmetric  behavior.  The  addition  of 
small  strakes  also  eliminates  the  asymmetry  quite  effectively. 

The  yawing  moment  generated  from  the  asymmetric  vortical  flow  pattern  was 
usually  large  enough  that  the  deflected  vertical  tail/rudder  could  not  produce  an  effective 
moment  to  compensate.  A number  of  forebody  vortex  control  techniques  have  been 
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reported  to  eliminate  the  lateral  instability  at  moderate  to  high  angles  of  attack.  Malcolm 
[27]  reviewed  several  experimental  efforts  which  include  movable  strakes,  blowing 
surface  jets,  blowing  and  suction  through  surface  slots,  suction  through  surface  holes  and 
miniatured  rotatable  tip  strakes.  The  idea  is  to  control  the  existing  forebody  vortex  so 
that  the  generated  yawing  moment  can  be  used  to  enhance  the  aircraft  maneuverability 
instead  of  degrading  it.  Gee,  Rizk  and  Schiff  [28]  numerically  investigated  the  effect  of 
forebody  tangential  slot  blowing  about  a F/A-18  aircraft  by  solving  the  thin-layer  Navier- 
Stokes  equation.  The  solutions  were  obtained  at  =0.243,  a =30.3°,  and 
/?e^=  11.0x10®.  The  original  Baldwin-Lomax  turbulence  model  was  used  through  the 
flowfield  except  the  forebody  region,  the  boundary-layer  diverter,  and  the  regions 
between  wing  flaps  and  fuselage.  In  the  forebody  region,  the  modified  Baldwin-Lomax 
model  by  Degani  and  Schiff  [15]  was  used.  In  the  boundary-layer  diverter  and  flap- 
fuselage  regions,  only  the  inner  eddy  viscosity  formulation  from  Baldwin-Lomax  model 
was  used.  The  predicted  vortex  breakdown  and  the  separation  lines  were  found  to  be 
located  consistently  with  the  flow  visualizations  obtained  in  the  wind-tunnel  tests.  The 
averaged  aircraft  yawing  moment  agreed  well  with  experimental  measurements. 


Research  Objectives 


The  ultimate  goal  of  this  study  is  to  investigate  the  complex  transonic  vortical 
flow  dynamics  through  a time-accurate  thin-layer  Navier-Stokes  computation.  The 
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specific  flow  problem  considered  is  a Mach  number  0.96  turbulent  flow  past  a Secant- 
Ogive-Cylinder  BoatTail  (SOCBT)  projectile  at  60°  angle  of  attack.  Due  to  the  lack  of 
experimental  measurements  on  the  geometrical  model  and  flow  conditions,  emphasis  is 
placed  on  the  fundamental  fluid  physics  rather  than  the  engineering-level  application,  i.e. , 
a qualitatively  rather  than  quantitatively  accurate  solution  is  sufficient.  It  is  well  known 
that  the  Reynolds-averaged  thin-layer  Navier-Stokes  equations  implemented  with  the 
Baldwin-Lomax  turbulence  model  for  ideal  gases  can  be  accurately  simulated  in  many 
high  speed  turbulent  flow  problems  as  long  as  the  flow  separation  is  not  strong. 
However,  in  a flow  with  large  crossflow  separation,  the  predicted  eddy  viscosity 
magnitude  had  been  found  to  be  too  high  to  induce  the  secondary  vortical  flow 
phenomena.  In  order  to  assure  that  the  mathematical  model  used  in  this  computation  is 
good  enough  to  represent  the  physics  of  the  vortical  flow,  the  numerical  simulation  with 
both  the  original  Baldwin-Lomax  turbulence  model  [16]  and  the  modified  version 
suggested  by  Degani  and  Schiff  [15]  is  conducted.  The  impact  of  the  strong  diffusion 
predicted  by  the  original  model  is  first  investigated.  Upon  obtaining  proper  Justification 
to  the  simulated  solution  with  the  modified  Baldwin-Lomax  turbulence  model,  a detailed 
description  of  the  formation  and  development  of  the  vortical  flow  is  addressed,  including 
a possible  onset  of  asymmetrical  vortex  shedding. 
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Figure  1.1  A schematic  vortical  flow  structure  in  a crossflow  plane. 


Figure  1.2  Simulated  surface  skin-ffiction  lines  on  an  unwrapped  projectile  surface. 


CHAPTER  II 

GOVERNING  EQUATIONS 


The  fundamental  equations  of  fluid  dynamics  are  based  on  the  universal  laws  of 
conservation  of  mass,  momentum  and  energy  [29].  Conservation  of  mass  results  in  an 
equation  called  the  continuity  equation.  Conservation  of  momentum  is  equivalent  to 
Newton’s  Second  Law  which  yields  a vector  equation  known  as  the  momentum  equation. 
Conservation  of  energy  is  identical  to  the  First  Law  of  Thermodynamics  giving  a fluid 
dynamics  equation  called  the  energy  equation.  In  addition  to  the  equations  developed 
from  these  universal  laws,  it  is  necessary  to  establish  additional  relationships  between 
fluid  properties  in  order  to  "close"  the  system  of  equations.  For  example,  the  equation 
of  state  relates  the  thermodynamic  variables  pressure  (p),  density  (p)  and  temperature  (7). 
For  a turbulent  flow  computation  based  on  the  Reynolds-averaged  equations,  a turbulence 
model  is  needed  to  relate  the  turbulent  flow  quantities  to  the  field  variables. 


Navier-Stokes  Equations  and  Equations  of  State 
Governing  equations  of  unsteady  compressible  fluid  dynamics  based  on  the  laws 
of  conservation,  or  the  Navier-Stokes  equations,  are  described  in  terms  of  the  Cartesian 
coordinates  as  follows.  All  body  forces  and  external  heat  addition  are  neglected. 
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Continuity  Equation:  Conservation  of  mass  results  in 

d{pu)  ^ bjpv)  ^ d(pw)  ^ Q (2.1) 

dt  dx  dy  dz 

where  p is  the  fluid  density  and  u,  v,  w are  the  velocity  components  in  the  Cartesian 
coordinate  of  the  x,  y,  z directions,  respectively. 

Momentum  Equations:  Conservation  of  momentum  in  the  x,  y and  z directions  can  be 
written  as 


d(pu)  d(pu^+  p)  dipuv)  d(puw) 


dt 
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(2.2b) 


d(pw)  + d(pwu)  ^ djpwv)  ^ d(pw^+  p)  = ^ ^ ^ (2.2c) 

dt  dx  dy  dz  dx  dy  dz 


where  p is  the  thermodynamic  pressure.  Ty  are  the  viscous  stress  tensor  components 
For  a Newtonian  fluid  they  have  the  form 
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dx  dy  dz  dx  i dx  dy  dz  ’ 


-\(du  dv  dWx  , ^ dv  2 /-dv  dw  d«, 

dx  dy  dz  dy  3 dy  dz  dx 


“ bx  by  dz  dz  3'^'  dz  dx  dy’ 
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in  which  Stokes  condition  [30]  that  X+2/3/x=0  has  been  assumed.  Here  jx  is  the 
coefficient  of  viscosity  (dynamic  viscosity)  and  X is  the  second  coefficient  of  viscosity. 
The  coefficient  /x  is  a function  of  temperature.  For  air  it  can  be  represented  accurately 
by  an  empirical  interpolation  formula  based  on  Sutherland  theory  of  viscosity  [31] 

T T + S 

M (2.3) 

where  the  subscript  oo  denotes  the  fluid  reference  state  and  the  constant  5’=198.6°R. 

Energy  Equation:  The  conservation  of  energy  yields  the  energy  equation 

^ ^ d[(E+  p)u]  ^ d[{E+  p)v]  ^ d[{E+  p)w]  ^ ^ ^ ^ ^ (2.4) 

dt  dx  dy  dz  dx  dy  dz 


with 
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= 


dT 


/3. 


dT 

UT  + VT  + WT  + K 

yy  yz 


dT 

|8,  = MT,^  + VT  , + WT,,  + K 

^ z zx  zy  zz 


Here,  E,  T and  k represent  the  total  energy  per  unit  volume,  temperature,  and  thermal 
conductivity  respectively.  The  heat  transfer  by  conduction  is  related  to  the  temperature 


gradient  through  the  Fourier  law  as 


q = - K VT 

The  thermal  conductivity  k can  be  related  to  the  viscosity  coefficient  (x  by  the  Prandtl 
number 


Pr  = !^Cp 

K 

where  Cp  is  the  constant  pressure  specific  heat. 


Governing  Equations  in  Vector  Form 

The  governing  equations  for  a compressible  fluid  flow  in  Cartesian  coordinates, 
Eq.(2.1),  (2.2)  and  (2.4),  can  be  written  in  the  vector  form  as 
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««♦  ^ (2.5) 

dt  dx  dy  dz  dx  dy  dz 

in  which  q is  the  unknown  vector,  E,  F,  G are  the  flux  vectors  and  E^,  and  are  the 
viscous  flux  vectors 
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K,  J 

(2.5a) 


Equations  of  State 

The  time-dependent  Navier-Stokes  equation  Eq.(2.5)  consists  of  a set  of  five 
equations  with  seven  unknowns  p,  u,  v,  w,  p,  E,  and  T.  Two  more  equations  which 
relate  p,  p,  T and  E can  be  obtained  from  the  equations  of  state  of  the  fluid.  For  perfect 
gas,  they  are 

p = pRT,  e = CJ  (2.6a) 

where  R is  the  universal  gas  constant,  e is  the  internal  energy  per  unit  mass  and  C„  is  the 
constant  volume  specific  heat.  Assuming  that  only  internal  energy  and  kinetic  energy  are 
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important,  the  total  energy  E can  be  written  as 


p[e  + 


+ w^)] 

2 


By  applying  the  following  relations, 


-^=7,  C^-C,=R  (2.6b) 

we  may  further  develop  the  two  equations  of  state  (Eq.  (2.6a))  into  the  following  forms 

p = (7-  \)[E-  -p{u^  + + w')]  (2.7a) 

T = Jzl[E-}.{u^+  w^)]  (2.7b) 

R p 2 

Now  the  variables  p and  T are  related  to  the  dependent  variables  p,  u,  v,  w,  E and 
therefore  the  Navier-Stokes  equation  can  be  solved  for  the  five  unknowns. 


Non-dimensional  Navier-Stokes  Equations 
The  governing  equations  are  put  into  non-dimensional  form.  The  advantage  of 
non-dimensionalization  is  to  reduce  the  number  of  parameters.  By  choosing  the  proper 
reference  quantities,  the  non-dimensional  variables  can  be  defined  with  an  asterisk  as 


27 


X 


* 


X 


t 


Lla^ 


u 


* 


w 


JL 

Moo 


♦ ^ ’T'  ♦ r~’  ♦ E # P * ^ ii 

K = — , T = —,  E = J,  p = — Ty  = 

°°  °°  Poo^oo  PooMqo  Poo^oo 

Here,  we  denote  the  free  stream  condition  by  the  subscript  oo . L is  the  characteristic 
length  scale  and  a^o  is  the  free  stream  speed  of  sound.  By  using  as  the  reference 
velocity,  the  Mach  number  does  not  appear  in  the  governing  equation,  thus  reducing  the 
number  of  characteristic  parameters.  The  reference  temperature  is  taken  as  the  free 
stream  temperature  which  can  be  obtained  by  the  relation  = yRT^,  i.e.,  = 

aJ/jR. 

By  substituting  the  non-dimensional  variables  into  Eq.(2.5),  we  have  the  non- 
dimensional  Navier-Stokes  equations  in  the  vector  form  as 


dq*  ^ dE*  ^ dF  * ^ 3G  * 
with  the  Reynolds  number  defined  as 


1 ^ dE;  dF;  dc;  ^ 

— ( — + — + — ) 
Re  dx*  dy*  dz* 


(2.8) 


Pco^ 

Moo 


The  vectors  q*,  E*,  F*,  G*,  E*,  F*  and  G/  now  have  forms  identical  to  those 
dimensional  vectors  in  Eq.(2.5a)  except  that  the  variables  here  are  non-dimensional  (with 
an  asterisk  *).  However,  the  explicit  expressions  for  the  non-dimensional  (3  are,  for 
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example, 


The  equations  of  state,  Eq.  (2.7a)  and  (2.7b),  are  also  given  in  the  non-dimensional 
forms 


In  subsequent  sections  the  asterisk  * is  left  out  from  the  non-dimensional  governing 
equations  for  convenience  unless  otherwise  noted. 


In  computational  fluid  dynamics  the  governing  equations  are  often  approximated 
by  a finite-difference  method.  Previously  we  have  expressed  the  governing  equations  in 
terms  of  a Cartesian  coordinate  system.  However  for  many  flow  problems  with  complex 
geometries  it  is  more  effective  to  use  a boundary-fitted  curvilinear  coordinate  system. 
By  a proper  coordinate  transformation,  the  boundary  surfaces  in  the  physical  plane  can 


p*  = (7- 1)[E*-Ip*(w*'+ v*'+ w’")] 


(2.9) 


T*  = y(y-l)[El-hu^^+  W*^)] 

p ^ 


(2.10) 


Transformed  Navier-Stokes  Equations 
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be  mapped  onto  rectangular  surfaces  in  the  transformed  domain.  Moreover,  the  physical 
grid  points  can  be  clustered  in  the  region  with  high  gradients,  yet  still  can  result  in  a 
rectangular,  uniformly  spaced  computational  domain.  This  is  especially  advantageous 
in  simplifying  the  finite-difference  equations,  although  a set  of  somewhat  more 
complicated,  transformed  Navier-Stokes  equations  is  obtained. 

Let  us  denote  {x,y,z,t)  as  the  physical  domain  and  introduce  a generalized 
curvilinear  coordinate  as  the  computational  domain  such  that 

^=|(x,y,z,f),  ri=  r)(x,y,z,t),  f(x,y,z,0,  r=  t 


where  t and  t represent  the  times  which  remain  unchanged  during  coordinates 
transformation.  By  the  chain-rule 


Jt 


A = A 

dx 

A = A 

dy 

A = A 

dz 


we  can  rewrite  the  non-dimensional  Navier-Stokes  equation  Eq.(2.8)  as 
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3t 


] 


] 


] 


] 


(2.11) 


The  metric  coefficients  etc.  appearing  in  the  transformed  governing  equations 


(2.11)  can  be  obtained  by  writing 


' dT  ' 

’ 1 0 0 o’ 

' dt  ' 

dx 
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% Vy  V, 

< 

> 

dy 

dt 
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1 

1 

dz 

(2.12a) 


Or,  in  the  other  way. 
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’ 1 0 0 o’ 
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_ Zt  Z^  Z^  Zj. 

y,  J 

Substituting  Eq.(2.12a)  into  Eq.(2.12b)  and  rearranging,  we  have 


’ 1 
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0 ’ 

’ 1 

0 0 

0 ' 

^x 
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Defining  the  Jacobian  of  the  transformation 


d(r,^,V,0 

d(t,x,y,z) 


10  0 0 


t f t t 


1 

det[A] 


Eq.(2.13)  gives  the  relations 


(2.12b) 


(2.13) 


fc  = - y,z^) , f,  = j(z^,  - z/P , f,  = - x^P  (2.1 4) 
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= -^/Ix-yr^y-Z^n 


^t  ^Sx  ySy  ^T^z 


with  the  following  identity  equations  of  metric  coefficients 


dr  J J dr}  J dlJ 


J dr)  J dfj 


A(^)+A(^)+A(iz)  = 0 

d^  J dr)  J d^  J 


ary  dr)  J di J 


where 


- 1 


Equation  (2.11)  now  can  be  expressed  as  [32] 


(2.15) 


(2.15a) 
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A(i)+  A[1(^<7+  ^e+^f+  ?g)i 

dr  J J ' ^ ^ " 

+ VE+  y]^F+  riP)] 

* ifi*  r/»  r,G)i 

= ?A)1 

+ —[i{r]  F + rjF  + r]G)] 
dp  J ^ ^ ^ 

+ A[l(fE  + tF  + tG  )]} 

in  which  the  metric  identities  Eq.(2.15)  have  been  employed.  It  can  be  further  written 
in  the  simple  conservation-law  form 


^ ^1  <216) 

3t  dr]  Re  d^  drj  d^ 


by  introducing  the  new  vectors 
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q = 1,  E = (^,q+ 


F - (r),q+  y)fi+  r]F+  r\fi)U,  G = {^fl+  ^p)U 


~ ~ (VxEy'^  Vpy)^J, 


G = (fE  + tF  + tG  )U 


where  the  vectors  q,  E,  F,  G,  Ey,  Fy  and  Gy  are  defined  by  Eq.(2.5a)  with  all  related 
quantities  being  non-dimensional.  Moreover,  if  we  define  the  contravariant  velocity 
components 


U = ^ = ^,+  tv  + t vv 


dt 


V = — = 71  + 7]  U+  71  V + 71  W 

IX-  V z 


dt 


W = ^^  - 


dt 


= r,+  r.«+  Lv+  Lw 


then  the  new  vectors  in  Eq.(2.16)  can  be  also  written  as 
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(2.17) 


Note  that  for  a fixed  grid  the  curvilinear  coordinates  rj  and  f are  invariant  with  the 
time  t.  It  follows  that  = y],  = — 0. 


Thin-Laver  Navier-Stokes  Equations 

For  high  Reynolds  number  flow,  the  diffusion  processes  along  the  body  surface 
are  often  negligibly  small  in  comparison  with  the  one  normal  to  the  surface.  Let  the 
body  surface  be  mapped  onto  a f=constant  plane  in  the  transformed  domain  then 

the  viscous  derivatives  as  well  as  the  thermal  diffusion  derivatives  in  ^ and  rj  directions 
are  neglected.  In  contrast  to  the  classical  boundary-layer  theory,  all  derivatives  in  the 
convective  terms  are  retained  to  make  the  simplified  equations  capable  of  dealing  with 
separated  and  reverse  flow  regions.  The  thin-layer  Navier-Stokes  equation  in  a vector 


form  is  written  as 
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where 


dr  9|  dr]  Re 
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Turbulence  Model 

A direct  simulation  of  a three-dimensional,  high  speed  (high  Re)  turbulent  flow 
is  practically  not  feasible.  The  main  problem  lies  in  that  the  time  and  length  scales  of 
the  turbulent  motion  are  so  small  that  the  number  of  grid  points  and  time  steps  required 
usually  exceeds  present  computer’s  capability  [29].  The  smallest  length  scale  known  to 
exist  in  a turbulence  motion  is  the  Kolmogorov  length  scale.  The  ratio  between 
Kolmogorov  length  scale  to  the  characteristic  length  scale  in  any  single  direction  for  a 
flow  problem  is  to  the  order  of  Re^'"^  [33].  In  three  dimensional  space  it  attains  a ratio 
of  cubic  order,  i.e.,  Re'^'‘^.  In  order  for  the  mesh  size  to  have  the  same  order  as  the 
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turbulence  length  scale,  a flow  problem  with  Re  = 10^  needs  a total  number  of  grid 
points  in  the  order  of  which  can  take  up  the  available  memories  in  today’s  most 
powerful  supercomputers.  Therefore,  in  a three-dimensional  turbulent  flow  simulation, 
we  need  to  increase  the  computing  resources  to  the  order  of  if  we  want  to  raise  the 
Reynolds  number  by  a factor  of  10  times,  not  to  mention  the  effect  of  choosing  a proper 
turbulent  time  scale. 

A feasible  alternative  is  to  solve  the  time-averaged  Navier-Stokes  equations,  also 
referred  as  the  Reynolds  averaged  equations  [29].  The  mass-weighted  Reynolds  averaged 
equation  in  this  study  takes  the  identical  forms  as  Eq.(2.5),  with  all  dependent  variables 
now  becoming  time-averaged.  For  laminar  flows  the  coefficient  of  viscosity  n=ni  is 
determined  by  the  Sutherland’s  law  (Eq.(2.3))  and  the  thermal  conductivity  k is  defined, 
assuming  a constant  Prandtl  number,  as  k,  = /Pr.  For  turbulent  flows,  the 
"effective"  viscosity  coefficient  ju  comprises  two  parts, 

where  ixi  is  the  laminar  viscosity  and  is  the  turbulent  (or  eddy)  viscosity  which  can  be 
represented  by  the  well-known  Baldwin-Lomax  turbulence  model  [16].  In  turbulent 
flows,  an  additional  transport  of  heat  is  caused  by  the  turbulent  motion.  Accordingly, 
the  thermal  conductivity  k has  the  form 

Cpix,  Cpfx^ 

K = Kj  K=  + 

' ' Pr  Pr^ 

in  which  Pr,  is  the  turbulent  Prandtl  number  with  the  value  of  0.9  for  air  [29].  The 
Baldwin-Lomax  model  is  a two-layer  algebraic  eddy  viscosity  model  which  relates 
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turbulent  eddy  viscosity  [x,  to  the  time-averaged  variables  in  such  a way  that  no  additional 


variables  are  introduced.  This  makes  the  Baldwin-Lomax  model  an  easy  one  to  apply 
and  costs  relatively  little  computer  time,  which  has  a great  advantage  in  three- 
dimensional  Navier-Stokes  computation.  As  a matter  of  fact,  the  Baldwin-Lomax 
turbulence  model  has  been  applied  extensively  in  many  computational  cases.  The  results 
are  accurate  and  well  justified  as  long  as  the  flow  separation  is  not  strong. 

The  Baldwin-Lomax  turbulence  model  was  proposed  after  Cebeci’s  two-layer 
model.  In  Baldwin-Lomax  model,  the  eddy  viscosity  ix,  is  given  by 


where  y is  the  normal  distance  from  the  solid  wall  and  y crossover  is  the  value  of  y which 
gives 

outer* 

In  the  inner  region,  the  PrandtI-Van  Driest  formulation  is  applied 


crossover 


crossover 


(2.19) 


(2.20) 


t ' inner 


where  / is  the  mixing  length  with  the  form 


I = ^y[l-exp(-  ^)] 

/\ 


CO  is  the  magnitude  of  the  vorticity  vector 


CO 


dy  dx  dz  dy  dx  dz 
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and  is  a law-of-wall  coordinate  defined  as 


(2.21) 


The  subscript  w designates  the  local  flow  properties  on  the  wall  and  is 

the  friction  velocity.  The  viscosity  coefficient  at  wall  consists  of  only  the  laminar 


viscosity  p,,  and  the  wall  shear  stress  \t„\  is  equal  to  k and  are  two 


constants  set  to  be  0.4  and  26,  respectively. 


The  Clauser  formulation  is  used  for  the  outer  region 


= KC,,f>F^F„^(y)  (2.22) 

where  K is  the  Clauser  constant  which  equals  to  0.0168.  = 1.6  is  an  additional 

constant,  and 
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y mn.x  I 


max  max 


wake 


= mm 


^wky  max  ^dif^^max 


The  quantities  and  are  determined  by  the  function 


F(y)  = y|co|[l-exp(- ^)]  (2.23) 

A 

where  F^  is  the  maximum  F along  a radial  grid  line  off  a solid  surface  and  is  the 
value  of  y where  F^  occurs.  = 0.25  and  represents  the  difference  between  the 
maximum  and  minimum  velocity  magnitude  along  the  radial  line.  Ff^^^iy)  is  the 
Klebanoff  intermittency  factor  given  by 
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pKletiy)  = [1+  5.5(-^f]-' 

y max 

where  = 0.3. 

In  application,  we  first  compute  the  vorticity  magnitude  | w | as  well  as  the  mixing 
length  / at  each  grid  point  along  a radial  line  from  the  surface.  The  computation  for 
ifi^inner  IS  thcn  Straightforward.  To  obtain  we  need  to  evaluate  which  is 

determined  by  the  profile  of  the  function  F(y).  The  final  eddy  viscosity  distribution  along 
this  radial  ray  is  obtained  by  taking  the  smaller  value  from  either  the  inner  or  the  outer 
eddy  viscosity  distribution. 

The  formulations  discussed  above  are  all  dimensional.  The  corresponding  non- 
dimensional  equations  for  Eq.(2.20),  (2.21)  and  (2.22)  are 

(f^)inner  = (2.20a) 


+ 


(2.21a) 


^p)outer  ^P^^cpP^wake^Kleb^^ 


(2.22a) 


in  which  the  vorticity  magnitude  |w|  and  F^^  are  normalized  by  aJL  and  a^L, 


respectively. 
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Modified  Baldwin-Lomax  Turbulence  Model 

As  discussed  in  Chapter  I,  the  Baldwin-Lomax  turbulence  model  is  not  suitable 
for  massive  crossflow  separations  in  the  high-a  flow  computation.  Degani  and  Schiff 
[15]  have  stated  that  when  the  crossflow  separation  is  relatively  small,  the  behavior  of 
the  function  F(y)  resembles  Case  (a)  in  Figure  2.1.  Only  one  well-defined  peak  occurs 
on  each  of  the  windward  and  leeward  sides  and  the  determination  of  and  is 
straightforward.  When  a is  larger,  the  characteristics  of  Fiy)  does  not  change  on  the 
windward  side.  However,  the  underlying  vortex  structure  above  the  leeward  side  can 
cause  a second,  larger  peak  as  shown  in  Case  (b).  A vorticity  magnitude  jump  takes 
place  when  the  radial  ray  crosses  the  vorticity  feeding  line  emanating  near  the  primary 
separation  point.  This  jump,  in  conjunction  with  a much  larger  distance  y from  the 
wall,  results  directly  in  the  second  larger  peak.  In  the  original  model  application,  the 
code  searches  along  the  outward  ray  and  determines  F^,  therefore,  picking  up  the 
second,  larger  peak  located  away  from  the  surface.  The  eddy  viscosity,  however, 
depends  only  on  local  flow  profiles  and  should  be  determined  within  the  attached 
boundary  layer  underlying  the  vortex  structures  on  the  leeward  side,  as  well  as  within 
the  windward-side  boundary  layer.  It  concludes  that  we  should  pick  up  the  first  peak 
within  the  attached  boundary  layer  and  exclude  the  second,  larger  one.  Indeed  the  choice 
of  the  larger  peak  can  result  in  a much  larger  computed  eddy  viscosity.  Consequently, 
the  turbulence  mixing  is  too  strong  and  the  details  of  the  flow  physics  and  phenomena 
in  the  thin-layer  can  be  distorted  or  washed  out. 

Degani  and  Schiff  then  suggested  the  following  modifications  on  the  Baldwin- 

Lomax  turbulence  model.  First  the  search  for  F^  is  performed  circumferentially  on  a 
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given  cross-section  from  the  windward  symmetric  plane  to  the  leeward  symmetric  plane. 
Along  each  radial  ray  the  code  searches  for  the  first  peak  it  encounters  then  cuts  off  the 
searching.  The  peak  is  considered  to  be  found  when  the  value  of  F(y)  at  the  next  radial 
grid  point  drops  more  than  10%  from  the  current  maximum  value.  Therefore,  choosing 
this  way  allows  us  to  obtain  a more  reasonable  computed  eddy  viscosity.  It  was  also 
mentioned  that,  when  the  location  of  the  vortex  feeding  sheet  appears  close  to  the  outer 
edge  of  the  attached  boundary  layer,  i.e.,  at  a very  short  distance  downstream  of  the 
crossflow  separation  point,  these  two  peaks  shown  in  F(y)  profile  tend  to  merge.  In  this 
case,  we  cannot  assume  a drop  of  more  than  10%  from  the  first  peak  before  reaching  the 
second  peak.  Accordingly,  a further  modification  was  made  in  such  a way  that,  on  each 
ray  a cut-off  distance  is  specified  in  terms  of  from  the  previous  circumferential  ray, 
i.e., 

= cy^(d-Ad)  (2.24) 

where  c is  a cut-off  constant,  and  6 increases  from  the  windward  symmetric  plane  to  the 
leeward  symmetric  plane.  When  a peak  does  not  appear  within  the  assigned  cut-off 
distance,  F^  and  y^  simply  take  their  values  from  the  previous  ray.  Determining  the 
constant  c depends  very  much  on  the  vortical  flow  structure  as  well  as  the  radial  grid 
discretization.  It  must  be  large  enough  for  the  code  to  detect  the  desired  peak,  but 
cannot  be  too  large  with  a possibility  to  pick  up  the  second  peak.  In  this  study,  the 
constant  c is  determined  to  be  5 through  numerical  experiments  which  are  described  in 
Chapter  V. 
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Figure  2.1  Behavior  of  the  function  F(y)  at  windward  (Case  (a)) 
and  leeward  (Case  (b))  locations. 


CHAPTER  III 
FLOW  PROBLEM 


In  this  chapter  the  flow  problem  and  geometrical  model  for  numerical  simulation 
are  defined.  The  governing  equations  are  the  thin-layer  Navier-Stokes  equations.  The 
domain  discretization  (grid  generation)  is  described  followed  by  a discussion  on  the 
imposed  boundary  conditions. 


Flow  Conditions  and  Geometrical  Model 
The  specific  flow  problem  considered  in  this  study  concerns  an  impulsively-started 
Mach  number  0.96  turbulent  flow  past  a Secant-Ogive-Cylinder-BoatTail  (SOCBT) 
projectile  with  sting  at  60-degree  angle  of  attack.  Defining  the  2.25-inch  diameter  on  the 
cylindrical  section  as  one  caliber,  the  real  projectile  model  consists  of  a 2.98-caliber 
forebody,  a 1.983-caliber  cylinder,  and  a one-caliber,  7°  boattail.  There  is  no 
measurements  available  for  assessing  the  numerical  simulation  of  the  flow  problem. 
However,  the  model  with  a 0.32-caliber  cylindrical  sting  attached  to  the  base  had  been 
investigated  extensively  by  the  U.S.  Army  Ballistic  Research  Laboratory.  The  surface 
pressure  coefficient  on  the  six-caliber  main  body  had  been  measured  for  the  angles  of 
attack  a at  0°,  4°,  6°  and  10°,  for  the  free  stream  Mach  numbers  at  0.91,  0.94, 
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0.96,  0.98,  1.10  and  1.20  [34].  The  Reynolds  number,  based  on  a free  stream  velocity 
V^=Q.96a^  and  the  diameter  of  the  projectile  cylinder,  is  760,000.  The  stagnation 
pressure  Po  = ^ atm  and  the  stagnation  temperature  Tg  = 580°R,  or  = 0.55  atm  and 
= 489. 6°R.  The  Prandtl  number  for  the  testing  condition  is  0.71  [35]. 

A simplified  projectile  model  for  numerical  simulation  is  shown  in  Figure  3.1. 
The  boattail  has  been  extended  for  another  1.82  calibers  onto  the  sting  to  avoid  the 
difficulties  in  a base  flow  simulation.  The  total  length  is  15.963  calibers.  Also  shown 
in  Figure  3.1  are  the  Cartesian  and  curvilinear  coordinates  of  the  flow  problem.  The 
Cartesian  coordinate  system  has  x as  the  projectile  longitudinal  axis,  y lies  in  the  right- 
hand-side  direction,  and  z in  an  upward  direction.  The  curvilinear  coordinates  rj  and 
f represent  the  streamwise,  circumferential  and  normal  directions,  respectively. 


Domain  Discretization 

The  three-dimensional  grid  system  employed  in  this  study  is  obtained  by  rotating 
a two-dimensional  planar  grid  with  respect  to  x-axis.  It  has  92  points  in  the 
circumferential  direction.  The  planar  grid,  which  has  110  points  in  the  streamwise 
direction  and  60  points  in  the  normal  direction  respectively,  is  a good  grid  for  the 
transonic  flow  problem  at  10-degree  angle  of  attack.  In  the  following  sections,  the 
integers  J,  K,  L are  used  to  denote  the  grid  points  in  rj  and  f directions,  respectively. 
JMAX,  KMAX  and  LMAX  represent  the  grid  dimensions  in  each  direction. 
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Planar  Grid  Generation 

The  two-dimensional  planar  grid  with  JMAX  = 1 10  in  ^ direction  and  LMAX=60 
in  f direction  is  generated  by  a hyperbolic  grid  generation  code  and  shown  in  Figure  3.2 
[36].  It  has  an  exponential  clustering  in  the  f direction  with  the  minimum  spacing  of 
2 X 10"^  beiwen  the  surface  (L=  1)  and  first  grid  lines  (L=2)  which  keeps  the  law-of-the- 
wall  coordinate  at  L=2  within  the  value  of  1 throughout  the  computation.  This 
spacing  is  indeed  small  enough  for  an  accurate  turbulent  boundary  layer  flow  simulation. 
The  outer  domain  is  approximately  24  calibers  away  from  the  projectile  and  is  large 
enough  for  an  accurate  Mach  number  0.96  turbulent  flow  simulation  for  a <10°  [37]. 
A detailed  planar  grid  distribution  near  the  surface  is  shown  in  Figure  3.3.  A few 
selected  J locations  and  their  local  diameters  are  identified  for  reference.  Note  that  the 
110  points  in  x-direction  are  distributed  such  that  95  points  are  located  on  the  main  body 
(nose  tip  to  boattail-sting  juncture).  Accordingly,  the  largest  grid  spacing  on  sting  is 
about  1.03  calibers.  It  is  realized  that  such  a grid  resolution  may  not  be  sufficient  to 
justify  the  computational  results  on  the  sting  segment.  Therefore  in  a subsequent  version 
of  the  planar  grid,  some  additional  21  points  are  added  to  refine  sting  resolution  for  more 
accurate  asymmetric  vortex  shedding  simulation.  A comparison  of  these  two  planar  grids 
on  the  sting  is  shown  in  Figure  3.4. 

Three-Dimensional  Grid  Generation 

The  three-dimensional  grid  is  obtained  by  rotating  the  planar  grid  with  respect  to 
the  projectile  axis  in  the  circumferential  rj  direction.  A cross-sectional  grid  topology  on 
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the  cylinder  is  shown  in  Figure  3.5.  Starting  from  6=0°  (K==l)  on  the  windward 
symmetric  plane,  a suitable  discretization  chosen  has  A0=1O°  for  O°<0<4O°,  A6=5° 
for  4O°<0<65°,  and  A6=3°  all  the  way  to  the  leeward  symmetric  plane  0=180° 
(K=47).  The  grid  is  symmetric  with  respect  to  the  angle-of-attack  plane  with  a total  of 
92  K-planes.  Previous  studies  show  that  A0=1O°  can  give  an  accurate  windward  flow 
simulation  while  A0=3°  in  the  crossflow  separation  region  is  sufficient  to  resolve  the 
secondary  vortical  flow  structure  on  the  leeward  side. 


Boundary  Conditions 

The  three-dimensional  grid  in  the  physical  domain  (x,y,z)  is  transformed  onto  a 
computational  domain  where  the  thin-layer  Navier-Stokes  algorithm  is  applied. 

The  transformed  computational  domain  with  a constant  discretization  is  shown  in  Figure 
3.6.  The  boundary  conditions  (i.e.,  the  solution  vector  [p,  pu,  pv,  pw,  Ef  ) on  the  six 
surface  planes  must  be  provided  for  the  computation  at  every  time  step  and  are  discussed 
as  follows. 

Solid  Surface  Boundary  Condition 

The  solid  surface  L=1  (or  f=0)  is  mapped  onto  the  plane  EFGH  on  which  the 
"no-slip"  condition  is  enforced  for  a viscous  flow  case,  i.e.,  u=v=w=0.  An  adiabatic 
wall  is  assumed  such  that  dT/dfl=0  at  L=1  where  ^=V(*/|  Vf|  is  the  unit  normal  vector 
on  the  surface.  The  density  p and  total  energy  E on  the  surface  are  found  by  first 
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considering  the  inviscid  normal  momentum  equation.  Knowing  that  the  contravariant 
velocity  components  U,  V and  W also  vanish  on  the  solid  wall,  the  three  inviseid 
momentum  equations  are  rewritten  from  Eq.(2.18)  as 

Now  multiplying  these  three  equations  by  f,,  ^y,  and  respeetively,  adding  them  up 
together  and  utilizing  the  metric  identities  Eq.(2.15),  we  have  the  resulting  equation  in 
the  form 


(IX  + l/X  (’J/X  ’l/X  ^zOPv'^  f/X  ^z^z^P^  ~ 0 

or 

^.Px^  ^yPy^  ^zPz  = ^ 

This  is  the  same  as  having  VfVp  — 0 or  dp/dfl=0  at  L=l.  By  ideal  gas  law  the  fluid 
density  is  solely  a function  of  pressure  and  temperature.  Therefore  we  also  have 
dp/dfl=0  at  L=1  for  an  adiabatic  wall  (dT/dH—O).  Accordingly,  p and  p at  L=1  can 
have  their  values  directly  from  L=2.  The  total  energy  E at  L=1  then  can  be  computed 
from  the  equation  of  state  Eq.(2.9). 
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Outer  Domain  Boundary  Condition 

In  Figure  3.6,  plane  ABCD  is  the  outer  domain  L=LMAX  on  which  the  free 
stream  conditions  are  assumed 

P = P^',  P = ^.v,cosa,  V = 0,  w = V sina,  E = E^,  p = p , T = T 

• ~ OO  ~ 00  ■ * 00*  Oo’j/^OO’  00 

where  a is  the  angle  of  attack. 

Circumferential  Boundary  Condition 

The  flow  domain  in  the  circumferential  direction  is  cut  open  and  "unwrapped" 
during  the  coordinate  transformation.  In  order  to  establish  a periodic  boundary  condition 
in  the  cut-off  plane,  two  extra  K planes  are  added  to  overlap  on  K=1  and  K=2  planes 
which  makes  KMAX=94.  The  boundary  conditions  at  K = 1 (plane  ABFE)  and 
K=KMAX  (plane  CDHG)  are  provided  with  the  solutions  at  K=KMAX-I  and  K=2 
from  a previous  time  step.  The  overlapping  is  chosen  at  the  windward  symmetric  plane 
in  order  to  minimize  the  adverse  effect  of  updating  periodic  boundary  conditions  in  the 
wake  region. 

Singular  Plane  J — 1 

J=  1 (or  I =0)  is  a straight  line  off  the  projectile  nose  in  the  physical  domain.  It 
is  mapped  onto  a singular  plane  ADHE  in  the  computational  domain.  Subject  to  the 
physical  meaning  that  the  solution  on  this  singular  plane  should  be  circumferentially 
invariant,  a zero-order  extrapolation  in  J direction  is  applied  to  obtain  the  individual  q\j=i 
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for  each  K and  L,  i.e., 


q\,=  ^ ilil  (3.2) 

■<17=1  3 

Equation  (3.3)  was  obtained  in  the  original  3-D  thin-layer  Navier-Stokes  code  by 
enforcing  a zero-slope  condition  at  J = 1.  It  is  understood  that  the  assumption  is  not  ideal 
for  an  angle-of-attack  flow  simulation.  The  local  accurary  can  be  improved  by  using  a 
linear  extrapolation 


y=i 


(3.3) 


where  Sj2  and  are  the  distance  between  the  points  at  J = 1 to  J=2,  and  J=2  to  J=3, 
respectively.  In  this  study,  however,  Eq.(3.3)  is  used  throughout  the  computation. 

A summation/averaging  procedure  then  is  taken 


KMAX-l 

E 9 


J=  1 


K=  2 


(3.4) 


j=  1 


KMAX - 2 


to  form  the  singular  plane  boundary  condition  at  J = I used  in  the  code. 


Downstream  Boundary  Condition 

J=JMAX  is  mapped  onto  plane  BCGF  in  Figure  3.6.  Due  to  the  fact  that  the 
sting  has  been  extended  quite  far  downstream  away  from  the  projectile  main  body,  the 
code  in  an  earlier  version  simply  sets  a zeroth-order  extrapolation  on  p,  pu,  pv  and  pw. 
For  example,  Pjmax  = Pjmax-i  and  etc.  The  total  energy  E then  is  computed  from 
Eq.(2.9),  assuming  that  the  pressure  is  equal  to  free  stream  pressure p„.  After  analyzing 
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the  intermediate  computational  results  near  the  end  of  projectile  sting,  however,  some 
unrealistic  pressure  and  temperature  contours  were  found  near  J=JMAX  by  enforcing 
the  free  stream  pressure  at  the  downstream  plane.  When  the  study  of  the  solution  on 
sting  segment  found  to  be  important  at  a later  time  for  r>5.0,  the  JMAX  = 131  grid 
simply  replaces  the  original  JMAX  = 110  grid,  along  with  a zeroth-order  extrapolation 
for  all  five  variables  p,  pu,  pv,  pw  and  E to  set  up  the  downstream  boundary  condition, 
i-e,  Psmax  = PsMAx-i  which  is  not  equal  to  anymore. 

Slow  Start 

The  flow  problem  is  subject  to  an  impulsively  started  flowfield.  To  avoid  the 
sudden  jump  between  the  initial  free  stream  velocity  at  L=2  grid  lines  and  the  viscous 
"no-slip"  condition  on  the  solid  boundary  L=1  lines,  a slow  start  technique  is  applied 
in  the  code  such  that  the  velocity  at  L=1  is  reduced  gradually  to  zero  over  the  first  30 
time  steps.  This  is  accomplished  by  setting  the  solution  vector  q on  the  solid  surface  as 

1.30  (3.5) 

with 

CO  = (10-150+  60')0^  e = (N-  l)/30 

where  N is  the  number  of  time  steps  over  the  first  30  steps  and  q^^  is  the  solid  surface 
boundary  conditions  previously  discussed.  The  function  co  has  a distribution  which 
increases  very  slowly  over  the  first  few  steps  and  recovers  to  one  smoothly. 
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SOCBT  Projectile  Model  with  Sting 


Figure  3.1  Projectile  model  and 

coordinate  systems. 


Figure  3.2  2-D  planar  grid  (JMAX  = 110). 
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Figure  3.3  Close-up  of  2-D  planar  grid  (JMAX  — 110). 
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Figure  3.4  JMAX  = 110  grid  (top)  and  JMAX  = 131  grid  (bottom)  on  the  sting. 
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T1(K)  6=0°.  K=1  (windward) 


Figure  3.5  Circumferential  grid  distribution. 
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CHAPTER  IV 

THIN-LAYER  NAVIER-STOKES  COMPUTATION 


In  this  chapter  the  three-dimensional  thin-layer  Navier-Stokes  TVD  code  used  in 
the  study  is  first  discussed.  Then  a number  of  algorithms  whieh  relate  the  solution  vector 
q = [p,  pu,  pv,  pw,  E\^  to  other  useful  flow  variables  such  as  the  pressure  and  shear 
stress  coefficients  are  described.  These  variables  are  used  to  present  the  results  in  the 
next  chapter. 

The  thin-layer  Navier-Stokes  Total  Variation  Diminishing  (TVD)  code  is 
developed  for  solving  the  transformed  thin-layer  Navier-Stokes  equation  (Eq.(2.18)). 
The  seheme  is  based  on  a non-iterative,  implicit  approximately,  factorized  finite- 
difference  algorithm  described  by  Beam  and  Warming  [38]  whieh  is  accurate  to  the  first 
order  in  time.  Instead  of  inverting  a very  large  system  of  equations  resulted  from  all 
three  spatial  direetions,  the  implicit  operator  is  factorized  into  three  one-dimensional 
processes.  The  TVD  code  is  a modification  of  an  original  version  of  ARC3D  thin-layer 
Navier-Stokes  code  developed  at  NASA  Ames  Research  Center  [39].  In  application  the 
TVD  code  can  be  viewed  as  the  same  as  ARC3D  code.  It  only  uses  a different  way  to 
eonstruet  its  numerical  dissipation  terms,  and  has  been  shown  to  be  more  efficient  and 
robust  for  the  transonic  turbulent  flow  simulation  [40,41]. 
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Both  schemes  had  been  extensively  described  in  various  applications.  In  the 
following  only  the  TVD  scheme  is  presented  without  further  elaboration.  For  simplicity, 
all  vectors  in  Eq.(2. 18)  are  dropping  their  hat  It  is  therefore  a good  reminder  that  the 
vector  q now  represents  [p,  p«,  pv,  pw,  E\^  and  Jq  is  the  physical  non-dimensional 
solution  vector. 


Thin-Laver  Navier-Stokes  TVD  Code 

The  transformed  unsteady  three-dimensional  thin-layer  Navier-Stokes  equation  ean 
be  approximated  by  the  following  finite-difference  equation  based  on  the  delta  form 
described  Beam  and  Warming 

L,L  L.Aq"  = R"  (4.1) 

where  the  difference  operators  L^,  and  and  the  right-hand-side  vector  /? " are 

= [/+  At5^A+  DjJ' 

L = [/+  At6  B+  D,r 
= {/+  Ar{b^C-Re-^b^{J-^MJ)]  + 

R"  = [-AT(b^E+  b^F+  b^G-Re-^b^S)+  D^f 

Here  Aq"  = q"~^^  - q"  is  the  delta-form  unknown  to  be  solved  and  / is  a 5x5  identity 
matrix.  A = dE/dq,  B = dF/dq,  C = dG/dq  and  M = dS/dq  are  the  Jacobian  matrices 
obtained  in  the  time  linearization  process  of  E,  F,  G and  S,  respectively  [42].  They  are 
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all  5x5  matrices  and  are  explicitly  expressed  in  Appendix.  5^,  and  6^  in  L^,  L and 
Lj.  are  second-order  central  difference  operators  which  keep  the  resulting  systems  of 
equations  to  be  solved  in  a block-tridiagonal  matrix  form.  In  /? " a fourth-order  central 
difference  approximation  is  used  for  the  convection  terms  {E,  F and  G)  in  order  to  keep 
the  convection  truncation  error  from  exceeding  the  one  of  the  viscous  terms.  For  the 
viscous  term  S the  second-order  central  difference  approximation  is  also  used.  Notice 
that  the  elements  of  matrix  M contain  derivatives  in  f direction.  For  those  terms  with 
the  form  appeared  in  b^J'^MJ  ),  they  are  differenced  as  [32] 


The  implicit  dissipation  operators  D,^,  D,^,  D,^  and  the  explicit  dissipation  term  appear 
in  Eq.(4. 1)  are  results  of  the  second-order  symmetric  TVD  scheme  [43,44].  These  terms 
are  artificially  added  in  order  to  damp  out  the  spurious  oscillations  that  occur  when 
dealing  with  flow  discontinuities.  They  are 


(4.2) 


(4.3) 


2 2 2 2 2 2 


2 2 2 2 2 2 
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where 


= diag[-max^(X“  ,)]Aj,2 

2 2 


ft.  „ 1.  = diag[- max^(X  i)]A„  i 


J,K+_,L 
2 


K + _ 
2 


K+_ 
2 


(4.4) 


^j,K,L.i  = diag[-max^(X“  _,)]A^^^i 

2 2 2 

Here  €j  = 1.5  and  e^=l  are  two  constants  which  control  the  amount  of  added  numerical 
dissipation.  Xj+,^”',  Xk+.^'”  and  Xl+>^"'  are  the  m"'  eigenvalues  (w=l,5)  of  the  Jacobian 
matrices  A,  B and  C,  respectively.  The  subscripts  J + '/i,  K + Vi  and  L+V2  imply  that 
the  Jacobian  matrices  are  evaluated  as  » Qj,k+'ax  and  respectively.  For 

example,  q}+^A,K,L  is  a symmetric  average  of  ^j+,,k,l  and 


T^(X)  is  an  entropy  correction  function  given  by 


(4.5) 


^(X) 


|X| 

(e'+  X')/2e 


,|X|  >e  = 0.25 
,|X|  <e  = 0.25 


and  diag(z)  denotes  a 5 x5  diagonal  matrix  with  the  diagonal  elements  equal  to  z.  The 
operators  , Ak+.^  and  Al+.,^  are  defined  in  the  form,  for  example. 


= u 


J+1,K,L 


U 


J,K,L 


The  Fs  in  Eq.(4.3)  with  the  subscript  J + V2,  K + Vi  and  L+Vi  are  the  eigenvector 
matrices  of  Jacobian  matrices  A,  B and  C,  respectively.  The  subscript  Vi  here  has  the 
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same  meaning  as  in  Eq.(4.5).  The  4»’s  in  Eq.(4.3)  are  5 X 1 vectors.  For  example,  the 
m'^  element  of  the  vector  is 


where 


(4.6) 


gj” 


S max[0 , min(  I a""  , 

^^2 


S = sign(a 

^ 2 


_ y-l  ^j  + i.k.l/i  + i.k.l  ^J.K.L'^J.K.L 

^^2  0-5(-/j.i,K,L+  •^J.K.l) 


y. 


The  T ’s  are  the  inverse  of  corresponding  eigenvector  matrices.  The  explicit  expression 
of  these  eigenvalues  and  eigenvector  matrices  can  be  found  in  Ref.  45. 


Solution  Method  and  Application 

The  thin-layer  Navier-Stokes  computation  is  a time-marching  process  which 
calculates  the  solution  vector  along  rj  and  f directions  at  a given  time  t from  the 
solution  ^ at  a previous  time  step.  The  computation  takes  17.2  Mwords  of  total  memory 
of  Cray  Y-MP  or  Cray  C-90  for  the  110x94x60  grid.  The  CPU  time  required  for 
each  100  time  steps  is  approximately  0.72  hours  on  Cray  Y-MP  and  0.37  hours  on  C-90, 
which  translates  into  4.2x10"^  and  2.14x10'^  seconds  per  time  step,  per  grid  point, 
respectively.  In  this  section,  a flow  chart  illustrating  the  time-marching  process  from  the 
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time  step  nio  n+\  is  first  presented.  Then  a few  algorithms  which  convert  the  solution 
vector  to  other  useful  flow  parameters  are  described.  Finally  the  evaluation  of  Jacobian 
J and  the  applicable  size  of  time  step  At  are  briefly  discussed. 

Algorithm  and  Flow  Chart 

Equation  (4. 1)  shows  the  delta  form  of  the  factorized  finite-difference  formula  for 
the  thin-layer  Navier-Stokes  code.  In  application  it  is  split  into  three  block-tridiagonal 
systems  in  each  spatial  direction  and  is  solved  in  the  following  sequence 

= R"  (4.7a) 

=Ag**  (4.7b) 

L^Aq"  = Aq^  (4.7c) 

qn*i  = q»+  Aq" 

where  Aq*  and  Aq**  are  two  intermediate  solutions.  A flow  chart  shown  in  Figure  4. 1 
presents  the  solution  algorithm  of  the  thin-layer  Navier-Stokes  code.  A brief  description 
of  some  subroutines  in  the  main  program  is  provided  as  follows. 

BC  :Set  up  the  boundary  conditions  described  in  Chapter  III,  i.e. , the  q vector 

at  J = 1 and  JMAX,  K = 1 and  KMAX,  L=1  and  LMAX. 

VISRHS  : Evaluate  the  viscous  right-hand-side  -Re  ^5^  S. 

: Evaluate  the  rest  of  the  right-hand-side  terms. 


RHS 
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FILTRX  :Fill  in  the  elements  of  the  block-tridiagonal  matrix  L^. 
BTRIX  :Solve  the  block-tridiagonal  system  Eq.(4.7a). 

FILTRY  :Fill  in  the  elements  of  the  block-tridiagonal  matrix  L^. 
BTRIY  :Solve  the  block-tridiagonal  system  Eq.(4.7b). 
VISMAT  : Evaluate  the  viscous  term  -Re^b^{J‘MJ)  in  L^. 

FILTRZ  :Eill  in  the  elements  of  the  block-tridiagonal  matrix  L^. 
BTRIZ  :Solve  the  block-tridiagonal  system  Eq.(4.7c). 


For  the  vortical  flow  simulation,  the  solution  process  is  divided  into  separated 
jobs  with  every  job  accounting  for  50  time  steps  computation  (i.e.,  NMAX=50  in  Figure 
4. 1).  In  each  job  the  thin-layer  Navier-Stokes  code  reads  in  the  grid  data  and  restart  data 
(solution  vector  q and  the  time  r at  the  end  of  previous  job)  at  the  head  of  the  main 
program.  To  the  end  it  writes  out  the  restart  output  ^-file  for  the  next  restart.  In  the 
50-time-step  computation,  the  code  computes  the  residual  at  each  10  time  steps  and 
records  them.  It  also  records  the  surface  pressure  coefficient  Cp  and  the  surface  shear 
stress  coefficients  Tj.^,  at  every  5 time  steps.  A time  history  with  a good  resolution 
is  especially  important  when  we  are  interested  in  the  transient  flow  development. 

The  Baldwin-Lomax  turbulence  model  is  employed  in  the  code  to  evaluate  the 
eddy  viscosity,  with  an  option  to  use  either  the  original  or  a modified  version  by  Degani 
and  Schiff.  In  the  modified  Baldwin-Lomax  turbulence  model  application,  it  is  necessary 
to  search  for  the  quantities  and  in  the  circumferential  direction  from  windward 
location.  To  eliminate  the  unwanted  asymmetry  in  the  current  study,  the  searching  is 
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made  to  start  from  the  windward  symmetric  plane  (0=0°  or  K = l)  and  end  at  the 
leeward  symmetric  plane  (0=180°  or  K=47)  through  both  sides  of  the  body. 

Residual 

The  numerical  algorithm  used  is  a time-marching  scheme.  The  steady-state 
solution  is  obtained  as  a time  asymptotic  solution  of  the  unsteady  Navier-Stokes 
equations.  In  this  case  we  have  Aq"  = 0,  or  ^?”  = 0 in  Eq.(4.1).  Therefore,  the 
convergence  criterion  to  the  steady-state  solution  can  be  represented  by  the  averaged  L2 
residual  of  the  vector  R"  over  the  internal  nodes  of  the  domain,  i.e., 

T 0 rpcirliiQl  = / J,K,L,M  2 V • / 

At  (JMAX  -2)(KM  ax  -2)(LM  AX  -2)  ^ 

where  J=2,  JMAX-1,  K=2,  KMAX-1,  L=2,  LMAX-1  and  M = l,  5.  In  general  the 
residual  drops  when  the  computed  solution  converges  to  a stable  flow  pattern.  It  is 
understood,  however,  the  steady-state  solution  may  not  exist  due  to  the  inherent 
unsteadiness  of  the  physical  flow,  for  example,  the  vortex  shedding  phenomena. 


Pressure  and  Shear  Stress  Coefficients 

By  the  equation  of  state,  we  may  have  the  pressure  p and  temperature  T as  well 
as  their  non-dimensional  quantities  p*  and  T*  expressed  as 
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w^)] 

(2.5) 

w^)] 

(2.6) 

w*^)] 

(2.9) 

^ + w*')] 

(2.10) 

The  pressure  coefficient  Cp  can  be  computed,  according  to  its  definition 


Cp  = 


P-P 
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p-p 
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where  p/p^a^^  is  the  non-dimensional  p*  shown  in  Eq.(2.9)  and  pjp^^a 


2 _ 


00 


p^RTJpJyRTJ  = 1/y.  Accordingly, 


2 yMl 


(4.9) 


To  compute  the  shear  stress  components,  for  example. 


. du  dv. 

to* 


with  its  non-dimensional  form  as 


.,3m*  ^ dv\  1 

Txy  = IJ-  (r— r 


ay  * ax  * 


where  = p^a„L/p„,  we  have  the  shear  stress  coefficient 
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On  solid  surface  the  non-dimensional  viscosity  coefficient  (i  consists  of  only  laminar 


viscosity  which  can  be  computed  by  the  Sutherland’s  law  (Eq.(2.3)) 
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where  T*  is  shown  in  Eq.(2.10)  and  >S'=198.6°R.  The  evaluation  of  velocity  gradients 
can  be  made  by  the  chain-rule,  for  example,  u^=u^^^+u^r}^+u^^^.  The  velocity  gradients 
M|,  and  Mj.  on  the  transformed  domain  are  computed  by  second-order  central  difference 
formula  at  the  interior  points,  and  forward/backward  difference  formula  at  boundary 
points. 


In  order  to  transform  the  stress  tensor  in  a Cartesian  system  to  the  one  in 
curvilinear  system,  we  have  [30] 


[r(|,»?,r)]  = [AY  [r{x,y,z)\  [A] 


(4.12) 


where  [t]  is  the  stress  tensor  and  [A]  is  a transformation  matrix  defined  as 


[T(|,t?,f)]  = 
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Each  individual  matrix  element  can  be  expressed  in  terms  of  the  metric  coefficients. 
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which  represents  the  directional  cosine  of  the  corresponding  axes.  For  example, 


i -e 
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Therefore,  the  shear  stress  components  Tj.^,  and  can  be  expressed  as 


where 


+ c^)  _ (ax^  + by^  + cz^ 
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It  should  be  noted  that  at  the  singular  plane  J==l,  all  derivatives  with  respect  to  the 
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circumferential  direction  rj  are  zero.  It  leads  to  that  Eq.(4.13b)  is  indeterminate. 
Therefore  the  shear  stresses  at  J = 1 are  not  computed  in  the  code. 

The  surface  points  with  zero  are  defined  as  the  crossflow  separation  points  in 
a crossflow  plane.  For  a given  circumferential  distribution,  the  zeros  can  be  located 
by  the  bi-section  method  together  with  second-order  Aitken  interpolation  and  are  used 
to  identify  the  crossflow  separation  and  attachment  points.  It  is  often  misleading  that  the 
crossflow  separation  points  are  where  the  flow  leaves  the  solid  surface,  i.e,  three- 
dimensional  flow  separation.  In  fact  the  three-dimensional  flow  separation  should  be 
associated  with  the  surface  singular  points  on  which  both  the  longitudinal  and 
circumferential  shear  stresses  become  zero. 


Aerodynamic  Force 

The  aerodynamic  force  vector  acting  on  the  projectile  can  be  obtained  by 
integrating  surface  pressure  and  shear  stresses  over  projectile  surface  consisted  of  panels. 
Except  for  those  panels  between  J = 1 and  2 which  are  triangles,  all  other  panels  are 
quadrilaterals.  Figure  4.2  shows  a typical  triangular  panel  and  a quadrilateral  panel. 
The  forces  due  to  pressure  p and  three  shear  stress  components  and  Zj-j.  on  each 

panel  can  be  assumed  to  be  the  geometrically  averaged  pressure  and  shear  stresses  acting 
on  the  panel  center  multiplied  by  the  panel  area  fl,  i.e.. 


T 


ft 


where  m—2>  for  triangular  panel  and  m=4  for  quadrilateral  panel.  With  the  ix,y,z) 
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coordinates  of  the  surrounding  nodal  points  given,  it  is  easy  to  compute  the  unit  vectors 
and  pointing  to  the  local  f and  rj  directions.  Also  n^=n^xn^  is  the  unit  normal 
vector  pointing  outward  from  the  panel  surface.  Then  the  resultant  force  vector  can  be 
expressed  as 


np 

F = = fJ  ^ fJ+  F^k 

1 = 1 

where  Qj  is  the  panel  area  and  np  is  the  total  number  of  surface  panels.  Finally  we  can 
decompose  the  force  components  to  get  the  lift,  drag  and  side  forces,  i.e.. 


lift  = FjCOs(a)  - F^sin(a) 
drag  = F^sin(o:)  + F^cos(a) 
side  force  = F 


Mach  Number 

Due  to  the  fact  that  all  velocity  components  are  normalized  by  free  stream  speed 
of  sound  a„,  it  is  easy  to  compute  the  Mach  number  based  on  as 

M = \/w^+  v^+ 

In  order  to  get  the  local  Mach  number  M based  on  the  local  speed  of  sound  a,  we  need 
to  correct  it  by  the  square  root  of  temperature  ratio,  i.e., 

M = M—  = M 
a 


where  T is  the  local  temperature. 
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Vorticitv 

Vorticity  of  a velocity  field  V is  defined  as 


where 


to  = VxV 
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The  magnitude  of  vorticity  then  is 
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Streamtraces  and  Skin-Friction  Lines 

It  is  often  useful  to  have  the  streamline  plot  to  study  the  flow  topology.  For 
unsteady  flow  problem,  there  exists  only  instantaneous  streamline,  or  the  streamtrace. 
Keep  in  mind  that  those  lines  are  different  from  the  steady-state  streamlines  or  particle 
traces. 

The  basic  idea  to  generate  a streamtrace  is  to  move  a meissless  particle  through 
a series  of  small  steps  in  the  direction  tangential  to  the  local  vector  field,  either  two  or 
three-dimensional  [46].  A predictor-corrector  integration  algorithm  is  used  for  the 
calculation  of  streamtraces.  At  each  intermediate  step  the  local  velocity  vector  can  be 
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obtained  by  a linear  interpolation  from  the  surrounding  nodal  vectors.  The  step  size  is 
controlled  which  typically  results  in  10  to  20  integration  steps  through  each  grid  cell, 
depending  on  the  grid  density  in  the  domain  interested.  The  starting  point  is  selected 
interactively  by  try-and-error  such  that  we  may  have  the  best  interpretation  of  the  flow 
structure.  By  reversing  the  velocity  vector  we  can  have  a backward  streamtrace. 

There  are  several  occasions  that  the  streamtraces  could  be  terminated.  The  first 
one  is  when  the  streamtrace  hits  the  preset  maximum  number  of  integration  steps.  We 
can  increase  the  number  of  maximum  steps  to  avoid  the  situation.  The  second  one  is 
when  the  streamtrace  goes  into  a stable  limiting  cycle  and  terminates  there,  no  matter 
how  large  the  number  of  integration  step  is.  The  third  one  is  when  the  streamtrace 
moves  outside  the  domain  boundary,  in  which  it  cannot  obtain  the  interpolated  velocity 
vector  and  is  terminated  on  the  boundary.  This  is  often  seen  near  the  solid  boundary  for 
a simulated  viscous  flowfield.  Theoretically  the  velocity  vector  very  next  to  a non- 
permeable  wall  should  be  tangent  to  the  local  surface.  In  a numerical  computation, 
however,  a slightly  inward  velocity  component  normal  to  the  surface  could  be  developed 
through  numerical  interpolation  or  integration  which  terminates  the  streamtrace  on  solid 
boundary. 

On  the  solid  surface  where  velocity  vanishes,  streamtraces  can  be  represented 
by  the  surface  skin-friction  lines  which  are  lines  tangent  to  the  two-dimensional  surface 
shear  stress  vector  field.  By  a proper  selection  of  starting  points,  the  surface  flow 
topology  can  be  clearly  presented. 
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Evaluation  of  the  Jacobian  J 

The  Jacobian  is  evaluated  according  to  Eq.(2.15a)  in  which  the  spatial  derivatives 
etc.  are  computed  by  second-order  central  difference  formula  at  interior  nodes  and 
forward/backward  difference  formula  on  boundaries.  It  is  noticed  that  the  Jacobian  at 
the  singular  plane  J = 1 is  computed  by  an  linear  extrapolation  from  the  Jacobian  at  J=2 
and  J=3. 

Integration  Time  Step  Ar 

In  the  time  marching  process,  a time  spacing  At  is  used  where  r is  the  non- 
dimensional  time  defined  in  Chapter  II.  If  we  consider  one  caliber  as  the  length  of  2.25 
inches  for  the  projectile  cylindrical  diameter  and  the  speed  of  sound  as  1085  (at  sea 
level  and  r^,=580°R),  one  unit  of  t is  approximately  equal  to  1.73  x 10^  seconds  in  real 
time.  To  reduce  the  cost  of  computation  resources,  it  is  desired  to  speed  up  the  time- 
marching process  as  fast  as  possible.  Although  the  linear  analysis  shows  that  the 
numerical  scheme  is  unconditionally  stable,  the  nonlinear  equations  are  subject  to  a time- 
step  limitation.  In  general  At  must  be  smaller  when  the  flowfield  is  subject  to  a rapid 
transition,  for  example,  higher  speed  or  higher  angle-of-attack  flow.  A most-often-seen 
error  suffered  from  a larger  At  is  that  the  computed  total  energy  E turns  out  to  be  less 
than  the  amount  of  dynamic  energy  such  that  the  internal  energy  e is  negative.  The 
temperature  then  is  negative  so  the  evaluation  of  the  laminar  viscosity  (Eq.(2.3)  or 
(4.11))  causes  a mathematical  error  (negative  square  root).  To  avoid  such  difficulties, 
the  maximum  At  allowed  is  in  the  order  of  10'^. 
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Figure  4. 1 Flow  chart  of  thin-layer  Navier-Stokes  code. 
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Figure  4.2  Quadrilateral  and  triangular  panels  used  to  computed  aerodynamic  forces. 


CHAPTER  V 

EFFECTS  OF  TURBULENCE  MODELS 


The  evolution  of  projectile  vortical  flow  based  on  the  original  Baldwin-Lomax 
turbulence  model  for  JMAX  = 110  grid  has  been  simulated  to  r=14.475  with  At=10'^ 
(Case  A).  The  flow  problem  has  also  been  solved  by  the  modified  Baldwin-Lomax 
model  with  JMAX=110  grid  for  t<6.8,  and  with  JMAX  = 131  grid  for  5.0<r<ll.l 
(Case  B).  Accordingly,  a tremendous  amount  of  three-dimensional,  time-dependent 
flowfield  data  has  been  generated  for  analysis.  It  is  understood  in  literature  that  the 
modified  model  can  give  acceptably  accurate  surface  pressure  distribution  for  subsonic 
and  supersonic  vortical  problems,  especially  in  the  secondary  flow  region  [15,47],  In 
this  chapter,  quantitative  implications  of  both  turbulence  models  are  explored  and 
discussed.  The  impact  of  strong  diffusive  process  of  the  original  model  upon  the  flow 
solution  on  the  solid  surface  as  well  as  in  the  wake  region  is  investigated.  The  effect  of 
cut-off  distance  to  the  computed  eddy  viscosity  distribution  in  the  modified  model  is  also 
examined  through  a series  of  numerical  experiments. 

It  should  be  mentioned  that  the  flow  pattern  remains  symmetric  with  respect  to 
the  angle-of-attack  plane  except  for  the  very  downstream  region  on  sting  with  jc>10. 
Therefore  only  half-of-the-domain  results  are  presented  throughout  most  of  the  chapter. 
When  referring  to  a longitudinal  J location,  the  JMAX  = 110  grid  in  Figure  3.3  is  used. 
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unless  otherwise  mentioned.  All  flow  quantities  presented  are  non-dimensionalized  as 
defined  in  Chapter  II  except  when  the  dimensional  units  are  explicitly  shown. 


Impacts  of  Excessive  Eddy  Viscosity 
On  the  Crossflow  Plane  J=54 

Results  simulated  at  t=9.325  for  the  J==54  cross-section  (near  the  middle  of  the 
cylinder  atx=4.042)  by  both  turbulence  models  are  first  selected  to  illustrate  the  effects 
of  excessive  eddy  viscosity  in  Case  A.  The  crossflow  vortex  pattern  is  found  to  be  in 
steady  state  for  both  cases.  Figure  5.1  shows  the  eddy  viscosity  contour  in  the  (y,z) 
plane  (rear  view  of  the  projectile).  Apparently  the  original  Baldwin-Lomax  model  has 
resulted  in  a much  larger  amount  of  eddy  viscosity  in  the  order  of  10^ ~ 10^,  compared 
to  the  near-zero  values  of  modified  model  solution  in  the  separated  region.  Figure  5.2 
shows  the  corresponding  instantaneous  streamtraces  in  a crossflow  plane.  A streamtrace 
is  a line  tangent  to  the  velocity  vector.  Here  only  the  Cartesian  velocity  components  v 
and  w are  used  to  generate  the  streamtraces.  Since  the  streamtraces  are  constructed  by 
using  small-step  integrations  over  the  velocity  field,  a small  change  of  transverse  velocity 
component  can  cause  a considerable  deviation  to  the  resulting  plot.  Taking  Figure  5.2 
for  example  where  the  v-component  in  y-direction  is  nearly  zero  on  the  leeward 
symmetric  plane  0=180°.  Even  the  starting  points  of  streamtraces  are  completely 
symmetric,  a small  disturbance  in  the  v-component  during  the  spatial  integration  has 
caused  the  streamtrace  to  move  significantly  to  one  side  of  the  symmetric  plane  as 
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shown.  It  should  be  emphasized  that  such  asymmetry  in  a streamtraces  plot  does  not 
necessarily  represent  a significant  quantitative  difference.  In  Figure  5.2  we  can  clearly 
see  that  the  primary  vortices  are  about  the  same  for  both  cases,  while  the  secondary 
separations  are  formed  close  to  the  leeward  side  surface  only  for  Case  B.  This  complies 
with  the  results  obtained  by  Degani  and  Schiff  [15]  that  the  excessive  computed  eddy 
viscosity  predicted  by  original  Baldwin-Lomax  model  can  suppress  the  onset  of  secondary 
separation.  Figure  5.3  and  5.4  show  the  vorticity  magnitude  and  F(y)  (Eq.(2.23)) 
contours.  Away  from  secondary  flow  regions,  the  excessive  eddy  viscosity  of  Case  A 
as  shown  in  Figure  5.1  seems  to  have  little  effect  upon  general  characteristics  of  the 
function  F(y)  as  well  as  vorticity  distributions  outside  the  thin-layer  region.  Figure  5.5 
and  5.6  show  the  Mach  number  and  Cp  contours.  Only  minor  shifts  of  the  contour  lines 
are  observed  except  near  the  secondary  separation  regions.  Although  the  flow  is  found 
supersonic  before  the  primary  separation  point  in  both  Mach  contour  plots,  the  sonic  line 
(M  = l)  close  to  surface  associated  with  the  decelerated  flow  from  M>1  to  M<1  is 
almost  parallel  to  the  local  flow  direction.  An  enlarged  portion  of  Mach  contour  for 
Case  B with  identified  thick  streamtraces  is  shown  in  Figure  5.5  for  illustration. 
Accordingly  the  velocity  component  normal  to  the  sonic  line  M = 1 is  subsonic  rather  than 
supersonic,  which  means  no  shock  wave  is  detected  in  the  computed  solution.  In  fact  the 
sonic  line  in  Figure  5.5  looks  very  similar  to  the  slip  line  observed  in  a supersonic 
boundary  layer  flow  which  sits  in  between  the  free  stream  supersonic  flow  and  the 
subsonic  boundary  layer  flow. 
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It  is  clear  from  Figure  5.1  to  5.6  that  the  use  of  a physical  (y,z)  plane  fails  to 
bring  out  quantitatively  the  secondary  flow  phenomena  close  to  the  projectile  surface. 
With  highly-stretched  grid  lines  near  the  body,  a transformed  (K,L)  plane,  on  the  other 
hand,  can  clearly  provide  important  flow  characteristics  in  the  thin  layer.  Figure  5.7  to 
5.12  show  the  corresponding  field  plots  in  the  (K,L)  plane,  compared  to  Figure  5.1 
through  5.6.  Only  one-half  of  the  cross-section  is  presented  since  the  flow  is  symmetric 
with  respect  to  the  angle-of-attack  plane.  Here  K = 1 and  K=47  represents,  respectively, 
the  circumferential  location  0=O°(windward  symmetric  plane)  and  0=1 80° (leeward 
symmetric  plane),  while  the  ordinate  L is  a node  along  an  exponentially  stretched  line 
from  the  surface.  The  normal  distance  at  node  L,  yL  measured  in  caliber  from  the  solid 
surface  (the  same  y used  in  the  Baldwin-Lomax  turbulence  model  described  in  Chapter 
II)  at  J=54  is  listed  in  Table  5.1  for  reference.  The  difference  of  eddy  viscosity  near 
the  surface  (e.g.,  L<30)  in  both  cases  is  now  clearly  presented  in  the  stretched  (K,L) 
plane  as  shown  in  Figure  5.7.  The  existence  of  secondary  flow  region  in  case  B is 
clearly  visible  in  Figure  5.8.  The  secondary  flow  separates  at  about  K=35  (0=144°) 
and  reattaches  at  K=28  (0=123°),  with  a size  about  0.183  calibers  in  span  and  0.009 
calibers  in  height  (L=23).  From  the  streamtraces  plot  we  can  easily  verify  the  crossflow 
topology  summation  rule  by  considering  the  entire  cross-section,  i.e.,  K = 1 to  93.  For 
example  in  Case  B of  Figure  5.8,  there  is  one  saddle  point  (S)  near  (K,L) =(47,46),  4 
nodal  points  (N)  at  two  primary  and  two  secondary  vortex  centers  (two  in  the  other 
corresponding  half  are  not  shown),  and  8 half-saddle  points  (S’)  corresponding  to  the 
crossflow  separation/attachment  locations  on  the  solid  surface  at  K = l,  18,  28,  35,  47, 
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Table  5.1  Normal  distance  node  L at  J=54 


L 

L 

Jl 

L 

Jl 

L 

1 

0 

16 

1.970x10' 

31 

4.995x10' 

46 

1.218x10" 

2 

2.000x10-' 

17 

2.457x10' 

32 

6.182x10' 

47 

1.505X10" 

3 

4.474  X 10-' 

18 

3.060x10' 

33 

7.650x10' 

48 

1.858x10" 

4 

7.536x10' 

19 

3.805x10' 

34 

9.467X10' 

49 

2.294x10" 

5 

1.132x10^ 

20 

4.729x10-' 

35 

1.172x10-' 

50 

2.830x10" 

6 

1.601x10-^ 

21 

5.870x10' 

36 

1.450x10' 

51 

3.490x10" 

7 

2.181x10-^ 

22 

7.282x10' 

37 

1.794x10' 

52 

4.303x10" 

8 

2.898x10-^ 

23 

9.030x10' 

38 

2.219X10-' 

53 

5.304x10" 

9 

3.786x10^ 

24 

1.119x10-' 

39 

2.746x10-' 

54 

6.535  X 10" 

10 

4.884x10"' 

25 

1.387X10-' 

40 

3.398x10' 

55 

8.048  X 10" 

11 

6.242X10-^ 

26 

1.718x10' 

41 

4.204x10' 

56 

9.905x10" 

12 

7.923x10-^ 

27 

2.127x10-' 

42 

5.201x10-' 

57 

1.218x10' 

13 

1. 000x10-' 

28 

2.634x10' 

43 

6.435  X 10' 

58 

1.496x10' 

14 

1.258x10' 

29 

3.260x10' 

44 

7.961  X 10' 

59 

1.837X10' 

15 

1.576x10' 

30 

4.036x10' 

45 

9.847x10' 

60 

2.402x10' 

59,  66  and  76  (K=59,  66  and  76  are  not  shown).  The  number  of  critical  points  satisfies 
the  summation  rule  Eq.(1.2)  such  that  4 + ‘/2(0)-l-'/2(8)=-l. 

In  Figure  5.9  we  notice  a vorticity  feeding  line  in  both  cases  that  emanates  from 
the  primary  separation  point  near  K = 18  and  ends  near  the  primary  vortex  center.  Unlike 
an  attached  boundary  layer  within  which  the  vorticity  magnitude  rapidly  drops  to  almost 
zero  near  the  edge,  here  the  vorticity  feeding  line  away  from  the  surface  has  created  a 
small  bump  of  vorticity  in  the  wake  region  along  L-direction.  This  bump,  combined  with 
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the  effect  of  a large  distance  y from  the  surface,  directly  results  in  a ridge  of  F(y)  along 
the  feeding  line  as  shown  in  Figure  5.10.  For  the  original  Baldwin- Lomax  turbulence 
model  application  the  code  picks  up  the  peak  and  corresponding  on  each  L-line, 
thus  producing  the  excessive  eddy  viscosity  distribution  shown  in  Case  A of  Figure  5.7. 
The  modified  model  in  Case  B chooses  the  first  peak  close  to  the  edge  of  boundary  layer 
and  results  in  a much  smaller  eddy  viscosity  magnitude,  especially  near  the  secondary 
flow  region  shown  in  Case  B of  Figure  5.7.  In  Figure  5.11  the  reduced  eddy  viscosity 
does  produce  shifts  in  some  Mach  contour  lines  in  the  secondary  flow  region.  Also  the 
Mach  number  near  primary  vortex  core  in  case  B seems  to  have  a smaller  value.  In 
Figure  5.12,  a slightly  broader  range  of  high  pressure  area  (Cp>-0J)  near  the  leeward 
side  surface  is  observed  from  K=23  to  32,  which  is  apparently  associated  with  the 
secondary  flow  separation  around  K=35.  In  general  the  large  eddy  viscosity  has  a 
negligible  effect  for  large  L (or  yj  outside  the  thin-layer  region. 

On  the  Surface  Flow  Properties 

In  addition  to  the  cross-sectional  view,  surface  flow  patterns  are  also  used  to 
illustrate  the  turbulence  models  effects.  Figure  5.13  and  5.14  show  the  surface  skin- 
friction  lines  and  corresponding  pressure/temperature  contours  at  r=9.325.  These  are 
side  views  of  the  projectile  model  covering  the  circumferential  domain  from  0°  to  180°. 
We  can  observe  distinctly  developed  primary  and  secondary  separation  lines  as  well  as 
a secondary  attachment  line  for  Case  B,  as  opposed  to  a much  simpler  skin-friction  lines 
pattern  in  Case  A.  Note  that  the  starting  points  for  skin-friction  lines  are  selected  by  try- 
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and-error  in  order  to  obtain  a proper  interpretation  to  the  surface  flow  topology.  It 
should  be  mentioned  in  passing  that  the  density  of  skin-friction  lines  shown  in  this  study 
does  not  possess  any  quantitative  meaning.  A further  investigation  of  the  solution  at  a 
later  time  shows  that  the  results  at  t=9.325  are  approaching  a steady-state  solution  upon 
most  of  the  surface  for  x<10,  except  for  a small  portion  near  the  cylinder-boattail 
juncture  in  the  leeward  separated  region. 

Again,  for  convenience  of  analysis,  skin-friction  lines  are  plotted  on  a stretched 
surface  (x,0)  which  is  obtained  by  unwrapping  the  circumferential  domain  to  a uniform 
span  for  all  longitudinal  locations.  Recall  that  the  ogive-cylinder,  cylinder-boattail,  and 
boattail-sting  junctures  are  located  atx=3.017,  5.0  and  7.82.  In  Figure  5.15  the  primary 
separation  line  sits  approximately  on  the  same  circumferential  location  in  both  cases, 
ranging  from  over  110  degrees  on  the  ogive  section  to  less  than  90  degrees  on  boattail. 
In  contrast  to  Case  B,  secondary  separation  and  attachment  lines  in  Case  A are  only 
observed  near  the  x>  10  region  on  sting.  In  Case  B we  also  see  a separated  region  in 
the  longitudinal  x-direction  in  addition  to  the  usually  observed  crossflow  separation  in  6- 
direction.  It  covers  approximately  from  x—4.S  (longitudinal  primary  separation)  to  5.5 
(longitudinal  primary  attachment)  and  from  0=90°  to  110°.  Details  of  the  longitudinal 
separated  flow  are  discussed  in  a later  section.  The  existence  of  such  an  axially  reverse 
flow  in  Case  B indicates  that  the  over-predicted  eddy  viscosity  in  Case  A does  produce 
a strong  turbulent  mixing  to  suppress  the  flow  separation,  although  a 7°  downslope 
boattail  is  indeed  prone  to  form  a longitudinally  reverse  flow.  For  Case  B in  Figure 
5.15,  a nodal  point  attracting  neighboring  skin-friction  lines  is  found  near  (x,0)=(4.8. 
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90°)  which  is  related  to  the  longitudinal  primary  separation  line  mentioned  above. 
Accordingly,  a saddle  point  connecting  the  longitudinal  primary  attachment  line  is  located 
near  (x,d)=(5.3,  80°)  to  balance  the  surface  flow  topology.  One  more  point  of  interest 
observed  in  Figure  5.15  is  the  off-centering  of  primary  attachment  line  (near  0=180°) 
after  x=8.5,  which  indicates  a loss  of  symmetry  with  respect  to  the  angle-of-attack  plane. 
This  can  be  noticed  in  both  cases,  implying  the  vortical  flow  simulation  implemented 
with  the  original  Baldwin-Lomax  turbulence  model  can  also  trigger  vortex  asymmetry 
regardless  of  the  excessive  eddy  viscosity  it  predicted. 

Figure  5.16  and  17  show  the  corresponding  surface  Cp  and  temperature  contours 
at  r=9.325,  with  a few  selected  streamtraces  added  for  reference.  On  the  windward  side 
before  the  primary  separation,  pressure  in  both  cases  drops  rapidly,  due  to  the 
accelerated  flow  in  a circumferential  direction  and  reaches  a minimum  value  right  before 
the  primary  separation  line.  In  the  separated  region,  a pressure  rise  is  observed  in  Case 
B on  seeing  the  secondary  attachment  line.  Situated  between  the  primary  attachment 
(near  0=180°)  and  secondary  attachment  lines  is  another  small  pressure  trough  that 
appears  right  before  the  secondary  separation.  Although  the  skin-friction  lines  pattern 
in  Case  A does  not  show  a secondary  flow  structure  over  the  projectile  mainbody, 
computed  surface  pressure  contour  turns  out  to  be  very  comparable  to  the  one  in  Case 
B.  The  surface  temperature  distribution,  on  the  other  hand,  is  more  correlated  to  the 
skin-friction  lines  pattern,  as  shown  in  Figure  5.17.  Primary  separation  lines  sit  right 
on  top  of  the  temperature  ridge  in  both  cases.  In  general  the  modified  Baldwin-Lomax 
turbulence  model  gives  a higher  surface  temperature  over  the  separated  region  by  a range 
of  20°  to  80°R,  especially  on  the  ogive  and  cylinder. 
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The  circumferential  surface  Cp  distributions  at  r=9.325  for  both  turbulence 
model  computations  are  presented  in  Figure  5.18  for  quantitative  comparisons.  Solid 
lines  represent  results  of  the  modified  model  (Case  B)  and  circles  are  for  the  original 
model  (Case  A).  Since  the  computed  solution  becomes  slightly  asymmetric  after  x=6, 
only  cross-sections  up  to  J=77  (x=  5.521)  are  selected  for  presentation.  Other  chosen 
J’s  are  J = 10  (x=0.554),  J=20  (jc=  1.555),  J=35  (jc=3.017),  J=54  (x=4.042)  and  J=68 
(x=5.0).  In  general  the  modified  Baldwin-Lomax  turbulence  model  only  yields  a slightly 
more  pressure  back  up  in  the  secondary  separation  region,  indicating  that  the  generated 
secondary  vortices  may  be  too  weak  to  make  significant  pressure  changes.  A larger 
deviation  of  circumferential  Cp  is  shown  at  J=68  and  77,  which  is  apparently  caused  by 
the  rather  complex  flow  pattern  in  the  separated  region  observed  in  Case  B (Figure  5.15). 

Effects  of  excessive  eddy  viscosity  in  Case  A are  also  observed  at  earlier  times 
during  the  vortical  flow  developing  stage.  Figure  5.19  to  21  show  the  surface  skin- 
friction  lines,  Cp  and  temperature  contours  at  r=2.0.  The  primary  separation  is  located 
further  to  the  leeward  side  compared  to  the  one  at  r=9.325.  The  secondary  flow  exists 
in  Case  A,  even  a tertiary  vortex  structure  is  found  on  boattail  in  Case  B.  Unlike  the 
solutions  at  t=9.325,  Cp  and  temperature  contours  at  r=2.0  show  large  disparity 
between  two  cases  due  to  very  different  surface  flow  patterns.  Figure  5.22  to  24  present 
the  same  field  plots  at  t=6.8.  The  flow  pattern  now  resembles  more  to  the  t=9.325 
case,  except  for  a difference  in  the  longitudinal  flow  separation  near  cylinder-boattail 
juncture  in  the  crossflow  separated  region.  Locally  the  flow  seems  to  be  still  in 
transition  where  only  a smaller  reverse  flow  area  in  jc-direction  is  observed. 


82 


Time  history  of  non-dimensional  forces  acting  on  the  projectile  model  is  computed 
and  presented.  Lift,  drag  and  side  forces  resulting  from  both  the  pressure  and  viscous 
shear  forces  are  shown  in  Figure  5.25(a).  Dash  and  solid  lines  represent  original 
Baldwin-Lomax  turbulence  model  (Case  A)  and  modified  model  (Case  B)  results, 
respectively.  The  tiny  circles  on  solid  lines  after  r=5.0  indicate  that  the  solution  is 
obtained  by  using  JMAX  = 131  grid  with  an  extrapolated  downstream  boundary  condition 
as  discussed  in  Chapter  IV.  Due  to  the  incorrect  simulation  near  the  end  of  sting  by 
previously  fixing  a downstream  boundary  pressure  in  Case  A as  well  as  in  Case  B for 
T<  6.8  with  JMAX  = 1 10  grid,  the  forces  presented  are  only  computed  up  to  half  of  sting 
at  x=11.88  (J=106  in  JMAX  = 110  grid,  or  J = 115  for  JMAX  = 131  grid).  In  Figure 
5.25(a)  only  little  difference  can  be  seen  regarding  the  lift  and  drag  forces,  especially 
when  they  converge  at  a larger  time.  Drag  forces  are  considerably  large  at  such  high 
speed  and  high  angle  of  attack.  The  asymmetric  vortex  shedding  on  sting  has  caused 
fluctuating  side  forces  as  shown  for  both  turbulence  models,  despite  of  certain  phase  shift 
is  observed.  The  oscillations  are  around  a zero  mean,  with  an  amplitudes  of  only  4% 
to  5 % of  the  lift  forces.  The  modified  model  seems  to  give  a smaller  non-dimensional 
period  of  1.9  rather  than  2.2  from  the  original  model,  or  in  terms  of  the  Strouhal  number 
as  0. 152  for  the  original  model  and  0. 176  for  the  modified  model.  The  Strouhal  number 
is  defined  as  fdlVo,,  where/,  d and  are  the  dimensional  frequency,  diameter  and  free 
stream  velocity,  respectively.  As  mentioned  in  Chapter  I the  Strouhal  number  for  an 
inclined  cylinder  flow  can  be  evaluated  as  0.21sin(a)  which  gives  the  value  of  0.182  for 


a=60°.  Notice  that  the  classical  Strouhal  number  of  0.21  is  based  on  the  measurements 
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performed  by  Roshko  [19,31]-  Apparently  the  modified  Baldwin-Lomax  model  has 
predicted  a more  comparable  shedding  frequency.  Figure  5.25(b)  shows  the  force  history 
due  to  viscous  skin-friction  only.  They  consist  of  barely  1 % of  total  forces  in  terms  of 
the  magnitude,  causing  little  change  to  the  overall  lift  and  drag  forces. 

On  the  Development  of  Crossflow  Separation 

As  mentioned  previously  in  Chapter  IV  that  the  crossflow  separation  or 
attachment  points  on  the  solid  surface  are  defined  with  their  circumferential  shear  stress 
component  as  zero.  In  Figure  5.26  the  time  history  of  crossflow  separation  and 
attachment  points  are  presented  at  six  selected  J cross-sections.  The  horizontal  axis  has 
a non-dimensional  time  t ranging  from  0 to  2 for  earlier  flow  development.  The  vertical 
axis  covers  a circumferential  domain  from  60°  to  300°,  with  the  absence  of  the 
windward  attachment  point  near  the  symmetric  plane  6=0°.  Hollow  circles  are  used  to 
represent  the  original  model  computation  Case  A,  while  solid  dotted  lines  stand  for  the 
modified  model  computation  Case  B.  Local  diameters  for  J = 10,  20,  54,  84,  95,  106  are 
0.261,  0.641,  1.0,  0.737,  0.32  and  0.32,  respectively.  Figure  5.26  shows  that  the 
primary  separation  initiates  at  the  leeward  symmetric  plane  0 = 180°.  The  primary 
separation  points  immediately  move  upstream  on  either  side  to  allow  the  primary  vortex 
recirculation  region  to  form  behind  the  body.  The  onset  of  primary  separation  is  found 
to  be  dependent  on  the  local  diameter.  For  a smaller  cross-section,  the  time  required  is 
less.  With  the  primary  vortices  further  developing  and  moving  away  from  the  surface, 
the  secondary  separations  and  attachments  are  initiated.  In  general  the  secondary 
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separation  region  in  Case  B covers  a larger  circumferential  span  than  it  does  in  Case  A, 
which  appropriately  reflects  the  effect  of  modified  Baldwin-Lomax  turbulence  model. 
The  secondary  flows  in  Case  B expand  in  size  as  time  grows;  accordingly,  pushing  the 
primary  separation  points  further  upstream  in  the  circumferential  direction.  For  a 
smaller  cross-section  near  the  sting,  development  of  tertiary  separations  in  Case  B is 
clearly  visible,  for  example,  0.8<t<1.2  at  J=95.  It  is  interesting  to  observe  a 
temporary  vanishing  of  secondary  flow  regions  in  Case  B at  J=95  and  106,  bringing  the 
primary  separation  points  slightly  back  toward  the  0=180°  plane.  The  secondary  flow 
regions  re-appear  after  a short  period  of  time  with  the  primary  separation  points  moving 
back  to  their  previous  upstream  locations.  Figure  5.27  shows  the  time  history  of 
crossflow  separation/attachment  development  we  have  simulated  for  the  above  J cross- 
sections.  The  modified  model  Case  B has  resulted  in  larger  and  more  sustained 
crossflow  secondary  flows  due  to  less  amount  of  eddy  viscosity.  At  J = 106  both  models 
predict  alternating  vortex  shedding,  regardless  of  the  much  larger  eddy  viscosity 
magnitude  in  Case  A.  The  vortex  shedding  phenomena  are  further  discussed  in  the  next 
chapter. 

In  Figure  5.28  the  onset  time  for  crossflow  separation  is  given  as  a function  of 
the  cross-sectional  diameter.  Since  different  cross-sections  on  the  ogive  and  boattail-sting 
may  have  identical  diameters,  solid  geometrical  symbols  and  thick  lines  are  used  to  high- 
light the  locations  on  boattail  and  sting.  Evidently  we  can  see  a linear  behavior  of  the 
onset  time  to  the  local  diameter.  All  symbols  lie  on  the  same  line  for  the  primary 
separation,  regardless  of  different  locations  on  projectile  or  different  turbulence  models. 
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As  a matter  of  fact  the  modified  model  does  not  play  a role  until  the  primary  separation 
is  formed,  therefore  the  onset  time  for  primary  separation  is  not  affected.  In  Figure  5.28 
the  modified  model  shows  an  earlier  onset  of  secondary  separation  on  the  ogive  and 
cylinder  (large  hollow  squares)  than  the  original  model  does  (large  hollow  triangles).  On 
boattail  and  sting  (small  solid  squares  and  triangles),  on  the  other  hand,  the  difference 
seems  to  be  limited,  which  implies  that  the  onset  of  secondary  separation  is  also  affected 
by  the  local  geometry  and  flow  characteristics  as  well  as  by  the  turbulence  models. 
Notice  that  some  of  original  model  results  for  the  onset  of  secondary  separation  (large 
hollow  and  small  solid  triangles)  are  isolated  with  dotted  lines  for  t near  1 and  up.  Due 
to  the  relatively  large  Ar  we  have  used  in  the  original  model  computation  within  that 
time  span,  the  results  inside  the  dotted-line  region  might  have  possessed  larger  truncation 
errors. 

With  the  help  of  time  history  plots  shown  in  Figure  5.26  and  27  we  can  easily 
visualize  the  impact  of  excessive  eddy  viscosity  on  the  crossflow  topology  at  any  instant 
time.  For  example.  Figure  5.29  and  30  show  the  streamtraces  by  two  different 
turbulence  models  at  J=54  cross-section.  In  Figure  5.29  streamtraces  are  plotted  at 
t=2.0  during  the  vortical  flow  development.  The  primary  separation  points  in  Case  B 
are  observed  at  slightly  upstream  locations  than  those  found  in  Case  A.  The  primary 
vortex  structure  is  larger,  with  secondary  separation  regions  clearly  visible  near  the 
surface.  In  Figure  5.30  the  crossflow  structure  becomes  nearly  steady  at  t=6.8,  in 
comparison  to  the  streamtraces  at  t= 9. 325  in  Figure  5.2.  Again  the  primary  vortices 
are  about  the  same  in  both  cases,  and  the  secondary  flows  can  be  only  seen  in  Case  B 
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due  to  the  less  amount  of  eddy  viscosity.  It  seems  that  the  effects  of  turbulence  models 
upon  the  primary  vortex  structure  (location  of  primary  separation  and  size  of  primary 
vortex)  are  more  significant  in  the  flow  developing  stage.  Figure  5.31  shows  the 
streamtraces  for  the  cross-section  J = 106  on  sting  at  an  earlier  time  (r=0.4)  before  the 
secondary  separation  occurs.  As  expected  the  primary  vortical  flow  structure  in  Case  B 
is  slightly  larger.  Figure  5.32  shows  the  streamtraces  at  t=0.9.  Unlike  in  Case  A 
where  simple  secondary  flow  regions  are  observed,  the  instantaneous  streamtraces  in 
Case  B exhibit  a more  complex  crossflow  structure.  Figure  5.33(a)  shows  the  enlarged 
secondary  flow  structure  for  Case  B at  J = 106  and  r=0.9.  Within  the  crossflow  plane 
there  are  6 nodal  points  (N),  3 saddle  points  (S),  8 half-saddle  points  (S')  but  no  half- 
nodes (N').  Note  that  a half-saddle  point  on  the  windward  symmetric  surface  is  not 
shown  in  either  Figure  5.32  or  5.33(a).  Also  in  Figure  5.33(a)  only  half  the  domain  is 
presented  due  to  symmetry.  A pair  of  nodes  and  a saddle  point  associated  with  the 
primary  vortex  structure  that  are  not  shown  in  the  enlarged  Figure  5.33(a)  can  be  located 
in  Figure  5.32.  Substituting  the  numbers  of  critical  points  into  Eq.(1.2)  we  obtain 
6-f  V2(0)-3-*/2(8)=-1  which  satisfies  the  crossflow  summation  rule.  Referring  to  Figure 
5.26  for  the  crossflow  separation  time  history,  however,  we  can  deduce  that  the 
crossflow  topology  should  consist  of  a tertiary  flow  region  inside  the  secondary  flow  near 
the  surface,  which  is  not  seen  in  Figure  5.33(a).  The  streamtraces  are  redrawn  in  one 
half  of  the  stretched  (K,L)  plane  and  the  tertiary  separation  region  is  clearly  identified 
in  Figure  5.33(b).  We  now  have  8 nodal  points,  7 saddle  points,  8 half-saddle  points  and 
4 half-nodal  points.  Again  critical  points  on  the  other  symmetric  half  and  the  half-saddle 
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at  K = 1 are  not  shown.  Although  the  number  of  critical  points  are  different,  the 
summation  rule  states  8 + V2(4)-7-V2(8)=-l  which  is  still  satisfied.  In  addition  to  the 
discrepancy  found  in  the  number  of  critical  points,  the  two  vortices  near  the  surface  in 
the  (y,z)  plane  Figure  5.33(a)  are  found  to  be  unstable  spiral  nodes.  In  the  (K,L)  plane 
of  Figure  5.33(b),  however,  they  are  apparently  stable.  It  is  mentioned  in  Chapter  IV 
that  the  streamtraces  are  obtained  through  a series  of  numerical  interpolation  and 
integration  in  the  vector  field.  An  equally-spaced  (K,L)  domain  in  general  is  believed 
to  yield  a more  accurate  streamtraces  distribution  than  a physical  (y,z)  domain  does  in 
a numerical  point  of  view  when  discrepancies  arise. 


Implication  of  the  Cut-off  Constant 

The  cut-off  constant  c in  Eq.  (2.24)  determines  an  appropriate  length  scale  across 
the  turbulent  boundary  layer  in  the  modified  Baldwin-Lomax  turbulence  model  to  reject 
excessive  eddy  viscosity  in  the  separated  region.  The  resulting  eddy  viscosity 
distribution  based  on  a given  cut-off  constant  is  found  to  be  very  dependent  on  the  local 
solution  profile  and  grid  characteristics.  Degani  and  Schiff  [15]  first  proposed  c=1.5 
for  their  supersonic  flow  case  (M  = I.8~3.0,  a=6.3°  -22.75°,  A0=2.5°  -5°)  to  match 
measured  pressure  distribution  and  surface  flow  direction.  After  that,  the  modification 
has  been  used  in  a number  of  applications,  yet  no  explicit  values  of  cut-off  constant  were 
reported  in  the  literature.  In  general  the  cut-off  distance  must  be  at  least  large  enough 
to  accommodate  the  boundary  layer  thickness  on  the  leeward  side.  If  the  cut-off  distance 
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is  too  small,  on  the  other  hand,  the  length  scale  is  limited  within  part  of  the  turbulent 
boundary  layer  and  the  resulting  eddy  viscosity  could  become  too  low.  As  a consequence 
the  leeward  suction  peak  generated  by  the  vortical  flow  structure  is  over-predicted  [47], 
In  this  section  the  author  wishes  to  examine  how  the  amount  of  cut-off  constant  can  affect 
the  resulting  eddy  viscosity  distribution.  The  computed  solutions  based  on  c=5  are 
provided  to  obtain  corresponding  eddy  viscosity  contours  for  various  cut-off  constants, 
according  to  the  modified  formulation  suggested  by  Degani  and  Schiff. 

The  effect  of  cut-off  constant  can  be  demonstrated  by  the  following  numerical 
experiment.  Based  on  the  simulated  flowfield  at  J=54  and  t=9.325  from  Case  B,  a 
number  of  constant  c’s  are  selected  to  compute  the  corresponding  eddy  viscosity 
coefficient.  Figure  5.34  shows  the  eddy  viscosity  contours  on  a partial  (K,L)  domain 
with  small  circles  representing  the  location  of  chosen  (i.e.,  y^).  Refer  to  Figure 
5.10  for  the  corresponding  F(y)  contour  plot.  Recall  that  if  the  peak  value  F^  is  not 
found  within  the  cut-off  distance,  values  of  F^  and  at  a given  K-plane  are  assumed 
to  be  the  same  from  the  previous  circumferential  K-location.  In  Figure  5.34  we  use 
shaded  circles  to  represent  such  locations.  The  primary  and  secondary  separations 
S2,  and  the  secondary  attachment  A2  are  also  shown  in  each  plot  on  the  horizontal  axis 
for  reference.  Streamtraces  pattern  for  this  case  can  be  referred  to  Figure  5.8.  It  is 
understood  that  flow  inside  the  thin  layer  has  a direction  from  attachment  points  A,  to 
separation  points  5,,  where  i=\  or  2.  In  Figure  5.34  the  predicted  eddy  viscosity 
magnitude  is  found  to  be  strongly  dependent  on  the  location  of  chosen  F^  and  y^.  In 
"upstream"  flow  regions  before  (1  <K<  18)  and  S2  (36<K<47),  the  computed  eddy 
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viscosity  contours  are  identical  for  different  cut-off  constants.  The  eddy  viscosity 
magnitude  "behind"  the  separation  points  S,  and  S2  grows  rapidly  as  the  constant  c 
increases.  The  terms  "upstream"  and  "behind"  used  here  are  relative  to  the  local  flow 
direction  in  the  thin  layer.  For  example,  a "upstream"  region  of  Sj  and  S2  are  located 
in  the  negative  and  positive  K-direction,  respectively.  When  c is  greater  than  6,  the  cut- 
off distance  simply  does  not  play  any  role  in  finding  the  peak  values  and  (no 
shaded  circles  are  present).  It  is  evident  that  a smaller  value  of  c (2.0  or  2.5)  is  effective 
in  reducing  large  amount  of  eddy  viscosity  between  Sj  and  S2.  When  c advances  from 
2.5  to  3.0,  the  predicted  eddy  viscosity  magnitude  "behind"  Sj  jumps  up  1 to  2 orders 
and  then  remains  constant  thereafter.  The  eddy  viscosity  magnitude  behind  S2,  on  the 
other  hand,  does  not  change  by  switching  c from  2.5  to  3.0,  but  picks  up  slowly  and 
steadily  with  increasing  c’s.  In  Figure  5.10  (the  F(y)  contour)  it  is  noticed  that  the  peak 
value  of  F(y)  for  Case  B is  elevated  off  the  surface  significantly  from  K = 19  to  24  behind 
primary  separation  Sj,  with  corresponding  L indices  as  22,  23,  26,  28,  29  and  31.  In 
Degani  and  Schiff’s  modification  the  cut-off  distance  on  current  K-plane  is  set  to  be  the 
constant  c multiplied  by  the  chosen  from  the  previous  K-plane.  Notice  that  F^  and 
y^  at  K=18  are  located  at  L=18.  For  a smaller  c such  as  2.5  the  cut-off  distance  at 
K = 19  has  a value  of  2.5xyjg  (y  at  L=18)  which  confines  the  search  of  F^  within 
L = 2 1 . Therefore  the  peak  F(y)  at  L = 22  is  not  reached  and  therefore  we  have  a shaded 
circle  at  K = 19.  Apparently  the  F^'s  at  K=20  to  24  all  fall  out  of  the  cut-off  distance 
that  assumes  the  value  from  previous  K in  a similar  situation  and  the  circles  are  shaded. 
With  a slight  increase  of  c to  3.0,  the  peak  F(y)  at  K = 19  is  reached  at  L=22.  The  cut- 
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off  distance  at  K=20  then  has  a value  of  3.0xy22  within  which  the  peak  F(y)  at  L=23 
can  be  easily  found.  The  effective  cut-off  distance  grows  rapidly  on  subsequent  K-planes 
due  to  an  increasing  ehosen  y^,  which  results  in  a jump  up  of  eddy  viscosity  distribution 
shown  in  Figure  5.34.  In  the  region  "behind"  secondary  separation,  however,  the  cut-off 
distance  is  limited  to  a relatively  small  value  from  previous  K-plane.  For  example, 
K=24  to  31  in  c=5.0  case.  In  fact,  the  shaded  circles  after  K=27  for  various  c’s  are 
all  caused  by  a local  peak  F(y)  with  the  approximate  value  of  0.14  around  K=25  to  27 
(Figure  5.10).  At  K=27,  has  a value  of  0.00473  loeated  at  L=20.  With  the  cut-off 
constant  c=5  the  seareh  of  F^  at  K=28  is  limited  within  yau-off—^  xO. 00473  =0.02365 
which  is  greater  than  y^j  but  smaller  than  y2s-  Although  in  Figure  5.10  the  first  peak  of 
F(y)  at  K=28  seems  to  be  located  at  L=26,  F(y)  at  L=27  did  not  drop  more  than  10% 
from  the  local  peak  value  at  L=26.  Therefore  y^  and  F^  are  taking  their  values  from 
K=27,  since  no  local  F(y)  peak  is  determined  for  L less  than  28.  Similar  conditions 
hold  for  K=29,  30,  31  and  32.  Unless  a cut-off  constant  c greater  than  5 is  specified, 
the  resulting  eddy  viscosity  magnitude  near  the  secondary  flow  region  is  restricted  to  a 
lower  value  as  shown.  It  seems  that  the  contour  plots  for  c=2.0  and  2.5  can  give  a 
smoother  eddy  viscosity  variation  along  the  circumferential  K-direction.  For  c=3.0  and 
up  the  sharp  drops  of  eddy  viscosity  magnitude  across  K=25  and  K=27~34  are 
apparently  inevitable  by  using  the  modified  Baldwin-Lomax  turbulenee  model  with 
Eq.(2.24).  Due  to  the  limited  computing  resources,  we  have  not  simulated  the  flow 
problem  with  different  c’s  in  the  modified  model  computation.  However,  it  is  believed 
that  by  effectively  reducing  the  eddy  viscosity  magnitude  in  the  separated  region,  a 


91 


qualitatively  correct  solution  can  be  accomplished  by  using  the  modified  Baldwin-Lomax 
turbulence  model  with  a cut-off  constant  c=5. 

Not  only  the  modified  model  can  give  quite  a different  /x, -distribution  at  different 
time  as  the  flow  develops,  the  dependency  of  computed  eddy  viscosity  to  the  cut-off 
constant  c is  also  varied  with  different  crossflow  topology.  Figure  5.35  shows  the  F(y) 
contour  plots  for  Case  B at  J=54  with  t=2.0,  and  J = 106  with  r=0.9.  The  top  graph 
represents  an  earlier  vortical  flow  structure  shown  in  Figure  5.29,  compared  to  the 
7=9.325  solution  at  J=54  in  Figure  5.2  and  5.8.  The  secondary  vortex  exists  with  its 
center  closer  to  the  solid  surface.  Accordingly,  the  first  peak  of  F(y)  in  the  secondary 
flow  region  near  K=30  is  located  at  L=22,  compared  to  L=27  in  the  r=9.325  case. 
The  bottom  F(y)  contour  in  Figure  5.35  represents  a flowfield  consisting  of  both 
secondary  and  tertiary  flows  (Figure  5.32  and  33).  The  peak  F(y)’s  associated  with 
secondary  and  tertiary  vortices  are  observed  near  L=25  and  L=  14,  respectively.  Figure 
5.36  shows  the  eddy  viscosity  contours  for  the  case  with  J=54  and  r=2.0  by  applying 
different  cut-off  constants  in  modified  Baldwin-Lomax  model,  based  the  solution 
computed  with  c=5.  Since  the  first  peak  F(y)  is  now  located  closer  to  the  surface,  a cut- 
off distance  with  c=5  is  large  enough  to  accommodate  the  first  peak  F(y)  in  the  separated 
region,  hence  no  cut-offs  (shaded  circles)  are  detected.  It  is  also  noticed  that  the  eddy 
viscosity  within  secondary  flow  region  jumps  up  more  than  one  order  in  magnitude 
around  L=20  when  c advances  from  4.5  to  5.0.  This  implies  that  eddy  viscosity 
magnitude  in  the  separated  region  is  very  sensitive  to  small  changes  of  c in  this  specific 
flow  condition.  For  the  other  case  at  r=0.9  and  J = 106,  however,  the  change  of  c from 
2.5  to  5.5  all  give  identical  eddy  viscosity  distribution  as  shown  in  Figure  5.37. 
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In  1993  Chaderjian  [47]  conducted  a study  for  the  comparison  of  two  Navier- 
Stokes  codes  in  simulating  high-incidence  vortical  flow.  The  specific  flow  problem 
solved  is  a Mach  number  0.2  flow  past  an  ogive-cylinder  with  a =20°  and  Rejj=5  x 10^. 
The  turbulence  model  used  was  based  on  the  Baldwin-Lomax  model  modified  by  Degani 
and  Schiff,  but  with  a different  way  to  specify  the  cut-off  distance  which  is  a constant 
throughout  the  circumferential  K-direction,  i.e.. 


cut-off 


- Cy 


max  I w 


(5.1) 


where  is  the  value  of  y^  at  the  windward  plane  of  symmetry  K = l.  C is  a user- 
specified  cut-off  constant  with  a suggested  value  of  3 to  5.  For  the  J =54  and  r=9.325 
case,  y^  at  the  windward  plane  of  symmetry  is  approximately  2.3  x 10'^  around  L=  17. 
A cut-off  distance  of  3 to  5 times  therefore  lies  between  L=22  to  26.  The  resulting  eddy 
viscosity  contour  is  similar  to  the  c=3  case  in  Figure  5.34,  but  with  a smoother  variation 
in  K-direction  and  smaller  computed  eddy  viscosity  near  and  S2.  Both  Eq.(2.24)  and 
Eq.(5.1)  used  in  the  modified  Baldwin-Lomax  turbulence  model  are  able  to  effectively 
reduce  the  amount  of  eddy  viscosity  in  the  separated  region.  Accordingly,  a realistic 
vortical  flowfield  is  obtained  without  sacrificing  the  advantage  of  using  the  simple 
algebraic  Baldwin-Lomax  model.  In  order  to  achieve  a quantitatively  accurate  solution, 
on  the  other  hand,  the  cut-off  constant  (either  c or  Q must  be  validated  through  carefully 


conducted  measurements. 
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Figure  5.1  Eddy  viscosity  contour  at  J=54,  r =9.325  - Case  A (left)  and  Case  B (right). 


Figure  5.2  Instantaneous  streamtraces  at  J=54,  r=9.325  - Case  A (left)  and  Case  B (right). 
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Figure  5.3  Total  vorticity  contour  at  J=54,  r=9.325  - Case  A (left)  and  Case  B (right). 


Figure  5.4  F(y)  contour  at  J=54,  t=9.325  - Case  A (left)  and  Case  B (right). 


Figure  5.5  Mach  contour  at  J=54,  t=9.325  - Case  A (left)  and  Case  B (right). 


Figure  5.6  Cp  contour  at  J=54,  r=9.325  - Case  A (left)  and  Case  B (right). 
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Figure  5.7  Eddy  viscosity  contour  at  J=54,  t=9.325  - Case  A (left)  and  Case  B (right). 
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Figure  5.9  Total  vorticity  contour  at  J=54,  r=9.325 
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Case  A (left)  and  Case  B (right). 
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Figure  5.10  F(y)  contour  at  J=54,  t=9.325  - Case  A (left)  and  Case  B (right). 
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Figure  5.11  Mach  contour  at  J= 54,  t =9. 325  - Case  A (left)  and  Case  B (right). 


Figure  5.12  Cp  contour  at  J=54,  t=9.325  - Case  A (left)  and  Case  B (right). 
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Figure  5.13  Surface  skin-friction  lines  (top),  Cp  (middle)  and  temperature  contours  (bottom) 

at  t=9.325  for  Case  A.  AC/>=0.1  and  Ar=20°R. 


Figure  5.14  Surface  skin-friction  lines  (top),  Cp  (middle)  and  temperature  contours  (bottom) 

at  7=9.325  for  Case  B.  AQ?=0.1  and  Ar=20°R. 
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Figure  5.15  Surface  skin-friction  lines  at  t=9.325  - Case  A (left)  and  Case  B (right). 
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Figure  5.16  Surface  Cp  contour  at  t=9.325  - Case  A (left)  and  Case  B (right). 
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Figure  5.17  Surface  temperature  contour  at  t=9.325  - Case  A (left)  and  Case  B (right). 


10 


Surface  Cp  Surface  Cp  Surface  Cp 


100 


Figure  5.18  Circumferential  surface  Cp  distribution  at  t=9.325 
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Figure  5.19  Surface  skin-friction  lines  at  r=2.0  - Case  A (left)  and  Case  B (right). 
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Figure  5.20  Surface  Cp  contour  at  t= 2.0  - Case  A (left)  and  Case  B (right). 
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Figure  5.21  Surface  temperature  contour  at  r= 2.0  - Case  A (left)  and  Case  B (right). 
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Figure  5.22  Surface  skin-friction  lines  at  r=6.8  - Case  A (left)  and  Case  B (right). 


<Xe  8 10  0 2 4X6  8 10 


Figure  5.23  Surface  Cp  contour  at  t=6.8  - Case  A (left)  and  Case  B (right). 
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Figure  5.24  Surface  temperature  contour  at  r=6.8  - Case  A (left)  and  Case  B (right), 
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Figure  5.25  (a)  Total  force  time  history,  (b)  Skin-friction  force  history. 
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Figure  5.26  Time  history  of  crossflow  separation/attachment  points  in  earlier  stage. 
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Figure  5.27  Time  history  of  crossflow  separation/attachment  points. 
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Figure  5.28  Onset  time  for  crossflow  separation  projectile  diameter. 
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Figure  5.29  Instantaneous  streamtraces  at  J=54,  r=2.0  - Case  A (left)  and  Case  B (right) 
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Figure  5.30  Instantaneous  streamtraces  at  J=54,  r=6.8  - Case  A (left)  and  Case  B (right) 
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Figure  5.31  Instantaneous  streamtraces  at  J = 106,  r=0.4  - Case  A (left)  and  Case  B (right). 


Figure  5.32  Instantaneous  streamtraces  at  J = 106,  t=0.9  - Case  A (left)  and  Case  B (right) 
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Figure  5.33  Instantaneous  streamtraces  at  J = 106  and  r=0.9  for  Case  B. 

(a)  In  (y,z)  plane  (upper),  (b)  In  (K.L)  plane  (bottom). 
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Figure  5.34  Computed  eddy  viscosity  contours  at  1=54  and  r=9.325. 
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Figure  5.35  F(y)  contours  for  J=54,  r=2.0  (top)  and  J = 106,  t=0.9  (bottom). 


112 


Figure  5.36  Computed  eddy  viscosity  contours  at  J=54  and  r=2.0. 
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Figure  5.37  Computed  eddy  viscosity  contours  at  J = 106  and  t=0.9. 


CHAPTER  VI 

TURBULENT  VORTICAL  FLOW  PHENOMENA 


In  the  previous  chapter  we  have  demonstrated  that  the  modified  Baldwin-Lomax 
turbulence  model  can  bring  out  physically  realistic  secondary  flow  phenomena  in  the  near 
wake  region.  Results  obtained  also  show  that  the  strong  diffusive  process  in  the  original 
Baldwin-Lomax  model  computation  has  little  qualitative  effect  upon  the  vortical  flow 
solution  in  the  outer  layer.  It  is  therefore  concluded  that  the  modified  model  can  indeed 
represent  turbulent  vortical  flow  physics  and  phenomena. 

In  this  chapter,  the  simulated  results  by  modified  Baldwin-Lomax  turbulence 
model  are  used  to  investigate  complex  three-dimensional  vortical  flow  phenomena.  The 
crossflow  streamtraces  at  selected  J cross-sections  are  first  presented  to  illustrate  the 
evolution  of  cross-sectional  vortical  flow.  The  surface  skin-friction  lines  and  related 
contours  are  shown  to  demonstrate  the  flow  development  on  projectile  body  up  tox=  10. 
A visualization  of  three-dimensional  flowfield  is  presented  with  the  aid  of  3-D 
streamtraces.  Finally,  the  vortex  shedding  phenomena  simulated  on  sting  section  are 

examined  through  a series  of  cross-sectional  streamtraces  and  surface  skin-friction  lines 
plots. 
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Cross-sectional  Vortical  Flow  Structure 

The  evolution  of  cross-sectional  vortical  flow  structure  of  the  impulsively-started 
60-degree  angle-of-attack  turbulent  flow  at  J=54  is  shown  in  a portion  of  (K,L)  plane 
in  Figure  6.1.  A saddle  point  which  separates  the  upward  and  downward  flow  at  K =47 
(0=180°)  is  observed  near  L=26  (y=0.0172)  before  the  cross-sectional  reverse  flow 
starts  at  t=0.36  (refer  to  Figure  5.26).  As  soon  as  the  crossflow  separation  point  begins 
to  appear,  a separation  line  originating  from  the  surface  is  formed  which  extends  further 
into  the  wake  region  at  a later  time.  The  flow  direction  on  the  separation  line  is  toward 
the  surface,  until  a stable  node  is  observed  at  r=0.45.  The  stable  node  moves  toward 
the  end  of  the  separation  line  at  r=0.60  and  is  to  develop  into  a stable  spiral.  Figure 
6.2  shows  an  enlargement  of  the  spiral  vortex  development  at  J =54.  The  spiral  remains 
stable  until  about  t=0.90  when  a stable  limiting  cycle  is  formed  within  which  the  spiral 
is  an  unstable  one.  As  the  unstable  spiral  is  fully  grown,  the  stable  limiting  cycle  simply 
disappears  and  the  unstable  vortex  pattern  takes  over.  It  is  understood  that  the  crossflow 
topology  can  be  transformed  from  a stable  spiral  to  an  unstable  one  through  the 
appearance  of  limiting  cycle,  and  vice  versa.  Visbal  and  Gordnier  [11]  suggested  that 
a stable  crossflow  spiral  vortex  is  associated  with  a positive  axial  velocity  gradient 
(vortex  stretching)  at  the  core,  if  the  flow  is  assumed  to  be  effectively  incompressible, 
while  an  unstable  vortex  is  present  due  to  the  negative  axial  velocity  gradient  (vortex 
compression).  Such  correlations,  however,  are  not  observed  in  the  current  flow  problem 
which  is  highly  compressible.  As  a matter  of  fact,  the  axial  velocity  gradient  du/d^ 
contour  plots  reveal  that  the  flow  is  decelerating  in  the  separated  region  at  all  times  ever 
since  the  stable  spiral  pattern  is  formed  at  r=0.7. 
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Depending  on  various  external  flow  conditions  the  evolution  of  primary  vortex 
can  be  very  different  from  the  one  at  J =54.  The  cross-section  J = 19  (x  = 1 .45)  on  ogive 
and  J=87  (x=6.59)  on  boattail  which  have  approximately  the  same  diameter  (D  =0.606 
and  0.608,  respectively)  are  selected  with  their  crossflow  streamtraces  presented  in 
Figure  6.3.  Unlike  the  crossflow  vortex  evolution  at  J =54,  the  one  at  J = 19  consists  of 
only  an  unstable  spiral  pattern  for  t>0.6.  The  velocity  component  at  leeward  symmetric 
plane  K=47  is  pointing  outward  from  the  surface  L=l,  due  to  a surface  curvature  from 
the  axially-increasing  ogive  diameter.  Immediately  after  the  formation  of  an  unstable 
vortex  at  r=0.6,  a pair  of  saddle  points  are  found  at  K=47  to  balance  the  crossflow 
topology.  The  existence  of  the  lower  saddle  point  connecting  the  upward  and  downward 
flows  due  to  surface  curvature  effect  was  also  reported  by  Ward  and  Katz  [48]  in  their 
experiments,  in  which  the  flow  visualization  of  an  inclined  body  of  revolution  was 
obtained  in  a towing  tank  with  a=40°~45°  and  /?eo=0.8xl0"-'2.7xl0^  The  vortex 
evolution  at  J=87,  on  the  other  hand,  looks  very  similar  to  the  one  at  J =54.  A limiting 
cycle  is  also  observed  at  t=0.6  which  transforms  the  stable  vortex  pattern  into  an 
unstable  one  at  a later  time. 

The  crossflow  streamtraces  for  t=6.8  and  1 1 . 1 at  J = 19,  54  and  87  are  presented 
in  Figure  6.4.  Despite  of  the  appearance  of  some  secondary  flow  regions,  the  primary 
vortex  remains  unstable.  Flow  patterns  at  J = 19  and  54  are  very  similar  from  t=6.8  to 
11.1,  indicating  that  a good  convergence  has  been  achieved.  The  flow  pattern  at  J=87, 
on  the  other  hand,  shows  a slower  convergence  due  to  its  rather  downstream  location  as 
well  as  the  inclined  7°  boattail.  The  secondary  flow  region  collapses  at  t=11.1. 
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although  it  is  still  seen  at  r=6.8.  Figure  6.4  also  shows  a growing  cross-sectional 
primary  vortex  in  the  longitudinal  J direction.  The  axial  velocity  component  («)  and 
crossflow  plane  velocity  magnitude  (square  root  of  v^-I-h^)  contour  plots  for  J = 19,  54 
and  87  at  r = 11.1  are  presented  in  Figure  6.5.  The  free  stream  axial  velocity  component 
has  a value  of  0.48  which  is  normalized  by  the  free  stream  speed  of  sound.  It  is  noticed 
that  the  axial  velocity  near  primary  vortex  core  decreases  dramatically  from  J = 19  to 
J=87.  In  the  region  above  primary  attachment  point  K=47,  on  the  other  hand,  the  axial 
velocity  slows  down  gradually  but  maintains  a relative  high  speed  compared  to  the  free 
stream  value  of  0.48.  The  slowing  down  of  crossflow  plane  velocity  magnitude  is  also 
observed,  especially  in  the  near  wake  region  under  primary  vortex  around  L=20.  Figure 
6.6  shows  the  corresponding  and  Cp  contour  plots.  In  contrast  to  the  mostly 

negative  du/d^  in  the  separated  region  at  J = 19  and  J=54,  the  one  at  J=87  is  found  to 
be  positive  near  the  primary  vortex  center  and  above  the  vanishing  secondary  flow 
region.  On  the  other  hand,  du/d^  is  apparently  negative  along  the  primary  separation 
line  at  J=87,  while  they  are  positive  at  J = 19  and  J=54.  The  Cp  contour  at  J = 19  is 
found  to  be  in  close  pattern  as  the  one  at  J=54,  only  with  a higher  surface  pressure  at 
the  windward  symmetric  plane  K = 1 and  a lower  one  in  the  wake  region.  The  surface 
pressure  at  J = 87,  however,  leaves  a rather  flat  distribution  in  the  wake  region.  Without 
a proper  circumferential  adverse  pressure,  the  development  of  secondary  separation  is 
simply  constrained. 

Figure  6.7  shows  a variation  of  crossflow  topology  along  the  projectile  at 
T=  1 1. 1.  The  size  of  primary  vortex  increases  in  axial  direction  on  ogive,  cylinder  and 
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boattail.  Near  the  boattail-sting  juncture,  the  primary  vortices  are  found  to  be  elevated 
away  from  the  surface  with  the  shrinking  cross-sectional  area,  while  small  recirculations 
are  formed  in  the  near  wake  region.  Eventually  the  primary  vortices  breakdown  in  the 
far  wake  region  and  the  ones  near  the  surface  simply  take  over  as  the  new  primary 
vortices.  The  crossflow  pattern  on  the  rear  sting  exhibits  alternating  vortex  shedding 
which  is  to  be  discussed  in  a separate  section. 


Development  of  Surface  Flow  Patterns 

The  development  of  surface  flow  pattern  is  illustrated  through  a series  of  skin- 
friction  lines  plots  as  shown  in  Figure  6.8.  Again  the  unwrapped  and  stretched  (x,6) 
surface  is  used  which  covers  the  projectile  from  x=0 ~ 10,  and  6—60°  ~ 180°.  The  flow 
at  r=0.05  remains  attached  on  most  of  the  projectile  surface  when  the  impulsive  flow 
is  initiated,  except  a small  portion  of  circumferentially  reverse  flow  is  observed  near  the 
nose  tip.  The  reverse  flow  soon  propagates  onto  the  ogive  along  the  longitudinal 
direction  at  t=0.2,  while  in  the  meantime  a primary  separation  line  is  formed  on  the 
sting.  Recall  that  the  onset  time  of  crossflow  primary  separation  was  found  to  be  linearly 
proportional  to  the  local  diameters  in  previous  discussion.  Eventually  the  primary 
separation  lines  on  ogive  and  sting  join  together  at  r=0.4,  and  the  secondary  flow  region 
starts  to  develop  on  the  ogive  section  as  well  as  on  the  sting  (t=0.8).  A tertiary 
separation  region  is  formed  at  r=0.8  right  behind  the  boattail-sting  juncture  (x=7.82) 
and  immediately  develops  onto  the  entire  sting  section  (t=0.9).  When  the  primary 
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separation  line  on  the  cylinder  moves  toward  the  windward  side,  the  tertiary  separation 
on  the  sting  section  begins  to  disappear  around  t=1.0.  The  tertiary  separation  located 
in  front  of  the  boattaill-sting  juncture,  however,  is  maintained  and  propagates  upslope 
along  the  boattail  at  a later  time.  The  secondary  separation  line  from  ogive  then  ties 
together  with  the  one  from  boattail  (t=  1.2)  and  the  secondary  flow  region  is  formed  on 
cylinder  section  (3  <x<5)  while  the  secondary  flow  on  sting  is  temporarily  gone  (t=  1. 1, 
1.2  and  1.4).  In  subsequent  development  the  secondary  attachment  line  promptly  moves 
upstream  in  the  circumferential  direction  together  with  the  primary  separation  line,  which 
allows  the  tertiary  flow  region  on  boattail  to  further  merge  into  the  secondary  flow.  The 
secondary  separation  and  attachment  lines  re-appear  on  the  sting  section  after  the 
secondary  flow  region  is  formed  on  the  cylinder  at  t=1.8.  It  is  noticed  that  the  initial 
flow  development  on  cylinder  is  mostly  under  circumferential  flow  action  rather  than  in 
the  longitudinal  direction.  This  can  be  seen  from  the  plots  before  t=1.2  on  which  the 
skin-friction  lines  are  joining  separation  lines  in  an  almost  perpendicular  direction.  The 
axial  flow  development  is  not  obvious  until  a jet  stream  is  formed  above  the  leeward 
plane  of  symmetry  after  the  secondary  flows  on  ogive  and  boattail  merging  together, 
which  creates  a more  complex  transient  pattern  before  a relatively  steady  surface  flow 
topology  is  formed  at  t=6.0. 

Figure  6.9  shows  the  development  of  surface  Cp  contours  in  a time  sequence 
associated  with  skin-friction  lines  plots  shown  in  Figure  6.8.  In  the  circumferentially 
upstream  region  prior  to  the  primary  separation  line,  the  pressure  declines  as  the  flow 
accelerates  toward  the  leeward  side.  Small  pressure  drops  are  found  longitudinally  on 
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both  the  ogive-cylinder  (x=3.017)  and  cylinder-boattail  junctures  (x=5.0)  due  to  local 
flow  expansions  at  the  corners.  On  the  other  hand  a pressure  rise  is  found  on  the 
boattail-sting  juncture  due  to  a flow  compression.  Pressure  troughs  are  located  slightly 
ahead  the  primary  and  secondary  separation  lines  in  the  circumferential  direction 
throughout  the  flow  development,  which  are  apparently  the  important  mechanisms  to 
accelerate  the  flow  before  reaching  the  separation  lines.  The  pressure  ridge,  on  the  other 
hand,  coincides  well  with  the  secondary  attachment  line  when  the  flow  pattern  is  getting 
steady  (t=4.8  and  up).  In  general  the  pressure  distribution  is  not  quite  sensitive  to  the 
apparently  weak  secondary  and  tertiary  flow  motions.  The  corresponding  surface 
temperature  contours  are  presented  in  Figure  6.10.  Similar  to  the  pressure  distribution, 
the  temperature  drops  from  the  windward  side  to  the  leeward  side  before  the  primary 
separation  line.  A temperature  ridge  is  observed  right  on  top  of  the  primary  separation 
line  for  r>3.0,  which  is  apparently  associated  with  the  pressure  rise  across  the 
separation  line.  Unlike  the  pressure  contour,  the  temperature  distribution  is  more 
sensitive  to  the  corresponding  flow  topology,  for  example,  the  tertiary  flow  development 
for  1.0<t<3.0. 

A number  of  transient  critical  points  observed  in  the  skin-friction  line  plots  are 
helpful  to  describe  the  formation  and  development  of  the  surface  flow  topology.  For 
example.  Figure  6.11  shows  the  earlier  flow  development  on  the  first  half  of  ogive 
section  for  1OO°<0<16O°,  in  terms  of  the  skin-friction  lines  and  surface  Cp  contour 
plots.  In  the  very  beginning  (r=0.05),  a saddle  point  is  found  close  to  the  ogive  tip  at 
135°  which  effectively  divides  the  surface  flow  into  four  regions.  To  the  left  and  top 
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of  the  saddle  point  are  the  longitudinally  and  circumferentially  reverse  flows, 
respectively.  A nodal  point  then  shows  up  at  r=0.1  bringing  a longitudinally 
downstream  flow  from  its  left.  The  pair  of  nodal-saddle  points  move  to  the  right  rapidly 
and  merge  together  at  r=0.4,  leaving  the  secondary  flow  developed  behind.  Eventually 
the  secondary  flow  is  formed  at  t-O.S  which  continues  to  propagate  onto  the  rear  part 
of  ogive  and  cylinder.  In  the  corresponding  pressure  contours  the  formation  of  two  low 
pressure  centers  associated  with  both  sides  of  primary  separation  line  and  a high  pressure 
center  associated  with  the  secondary  attachment  line  are  clearly  observed  from  r=0.2  to 
0.4.  The  high  pressure  center  cuts  into  those  two  low  pressure  regions  at  a later  time 
which  sets  up  an  overall  circumferential  pressure  distribution  we  have  seen  in  Figure 
5.18. 

Another  example  is  the  evolution  of  longitudinally  reverse  flow  shown  in  Figure 
6.12.  When  the  general  flow  pattern  is  nearly  stable  on  the  projectile  surface,  a 
secondary  flow  region  in  the  longitudinal  direction  near  the  cylinder-boattail  juncture 
starts  to  develop.  The  transformed  domain  presented  has  the  ranges  4<x<6  and 
60°  <6<  180°.  Recall  that  at  x=5  the  projectile  begins  to  see  a 7°  downslope  boattail 
from  the  cylinder  section,  which  is  in  favor  of  forming  a longitudinal  recirculation 
pattern.  Indeed  the  reverse  flow  is  observed  at  t=7.0  between  the  circumferential 
primary  and  secondary  separation  lines.  After  going  through  a rather  complex  transition 
only  the  reverse  flow  between  primary  separation  line  and  secondary  attachment  line  is 
maintained  at  t=  1 1.0.  The  remaining  two  pairs  of  nodal-saddle  points  formation  allows 
the  recirculation  area  to  reside,  while  in  the  meantime,  sets  up  a proper  link  between  the 
crossflow  primary  separation  and  secondary  attachment  lines. 
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The  surface  flow  pattern  appears  to  be  quite  steady  at  the  end  of  the  modified 
Baldwin-Lomax  model  computation  (t=  1 1. 1).  Figure  6.13  shows  the  skin-friction  lines 
and  the  Qj/temperature  contours  at  t=10.6  and  11.1.  Very  much  resemblance  is 
observed  except  on  the  sting  section  where  the  flow  pattern  is  in  transition  associated 
with  the  alternating  vortex  shedding.  Compare  Figure  6.13  to  the  skin-friction  lines  at 
t=6.8  (Figure  5.22)  and  t=9.325  (Figure  5.15),  we  can  see  the  diminishing  of 
crossflow  secondary  flow  structure  for  5.5  <x<  10.  In  fact  at  r = 1 1 . 1 there  is  no  cross- 
sectional  secondary  flow  region  detected  in  the  range  6.78<x<9.05  (J=88  to  104). 


Three-Dimensional  Streamtraces 

The  three-dimensional  instantaneous  streamtraces  are  presented  in  a time  sequence 
to  visualize  the  vortical  flow  development.  Due  to  the  hardware  limitation,  the  plotting 
is  made  by  using  only  upper  half  of  the  domain  with  every  other  grid  points  included  in 
all  three  directions.  In  Figure  6.14  the  particles  are  released  at  J=4  (x=0.086)  right 
above  the  surface  (L=2),  from  13  circumferential  locations  ranging  from  0=90°  to  270° 
with  A0=  15°.  At  earlier  times  when  the  axial  flow  has  not  well  developed,  for  example, 
T=  1.0  and  2.0,  the  streamtraces  are  only  able  to  reach  part  of  the  projectile  then  elevated 
from  the  surface  in  a direction  parallel  to  the  free  stream  flow.  When  the  streamtraces 
arrive  onto  the  boattail  at  a later  time,  they  are  first  found  to  stay  close  to  the  surface 
then  take  off.  Referring  back  to  Figure  6.9  for  the  pressure  distribution  from  t=3.0  to 
6.0  and  Figure  6.13  with  r=10.6  and  11.1  reveals  that  the  surface  pressure  on  boattail 
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exhibits  a low  value  of  -0.5  at  t=3.0  then  gradually  goes  up  to  -0.1  at  t=11.1.  The 
low  pressure  suction  formed  during  the  initial  boattail  flow  development  around  r=3.0 
is  believed  to  attract  the  streamtraces  close  to  the  surface  as  shown.  Figure  6. 15  shows 
the  streamtraces  at  t=11.1,  with  the  particles  released  from  J=35,  54,  68  and  77, 
respectively.  The  numbered  labels  1 to  6 attached  to  the  individual  streamtrace  imply 
that  the  particles  are  released  at  0=90°,  120°,  150°,  210°,  240°  and  270°,  respectively. 
Note  that  the  projectile  body  is  made  transparent  in  order  to  see  the  streamtraces  on  the 
other  side.  From  these  plots  it  is  evident  that  the  particles  released  near  0=90°  (negative 
y-direction)  and  270°  (positive  y-direction)  have  carried  mostly  circumferential 
momentum,  therefore  depart  from  the  surface  immediately.  The  particles  released  near 
the  leeward  symmetric  plane  (number  3 and  4),  on  the  other  hand,  have  more  axial 
momentum  to  transport  the  particles  downstream.  At  the  cross-section  J=77  the  traces 
number  1 and  6 are  missing  because  of  the  spatial  integration  based  on  local  velocity 
vectors,  which  has  caused  the  trajectories  to  fall  outside  of  the  current  plotting  domain, 
therefore  terminates  the  tracing. 


Vortex  Shedding  Phenomena 

As  mentioned  in  previous  discussion  that  the  simulated  vortical  flow  pattern  lost 
the  symmetry  with  respect  to  the  angle-of-attack  plane  and  the  alternating  vortex  shedding 
is  observed  on  the  sting  section.  Although  the  mechanisms  that  lead  to  the  vortex 
shedding  phenomena  are  not  well  understood  at  the  present  time,  a number  of  suggestions 
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have  been  made  by  Degani  [22]  to  relate  the  asymmetry  to  instabilities  in  the  flowfield. 
In  order  to  simulate  the  von  Karman  vortex  shedding  street  observed  behind  a 2-D 
cylinder  flow,  a symmetry-breaking  perturbation  must  be  applied  for  a short  period  of 
time  to  break  the  symmetric  solution.  The  perturbation  then  is  removed  and  the  flow 
was  able  to  advance  to  a self-sustained  oscillatory  solution.  This  is  called  an  absolute 
instability  of  the  symmetric  flow,  which  means  the  disturbance  will  grow  in  time  and 
space  once  it  is  introduced.  Similarly  in  a 3-D  ogive-cylinder  flow  simulation  performed 
by  Degani  at  very  high  angles  of  attack  (60°  and  80°),  a small  temporal  disturbance  of 
finite  duration  was  found  sufficient  to  trigger  the  unsteadiness  that  evolved  to  an 
alternating  vortex  shedding  in  the  wake  region  of  the  cylindrical  afterbody.  The  vortical 
flow  was  instantaneously  highly  asymmetric,  but  was  fluctuating  around  a symmetric 
mean,  revealed  by  the  total  side  force  time  history. 

In  this  study  it  is  believed  that  the  asymmetric  vortex  shedding  observed  is  a 
result  of  absolute  instability.  Due  to  the  one-sided  sweeping  from  K=I  to  KMAX  in 
solving  block-tridiagonal  system  along  rj-direction  (Eq.(4.7b))  in  the  thin-layer  Navier- 
Stokes  code,  an  asymmetric  truncation  error  is  introduced  in  the  transient  solution  which 
triggers  the  instability.  The  Strouhal  number  of  the  simulated  vortex  shedding  is  found 
to  be  quite  comparable  to  an  empirical  value  0.182  (as  discussed  in  previous  chapter) 
which  gives  a good  quantitative  assessment  to  the  numerical  accuracy. 

The  onset  of  vortex  asymmetry  is  first  observed  on  the  surface  skin-friction  line 
pattern,  which  is  especially  visible  near  the  leeward  symmetric  plane.  In  Figure  6. 16  the 
primary  attachment  line  is  shown  in  a time  sequence  on  a transformed  (x,6)  surface  with 
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160°  <0  <200°  for  the  overall  projectile  length.  The  oscillation  pattern  on  sting  section 
appears  like  the  undulating  of  a string  from  one  constrained  end  (boattail-sting  juncture 
at  x=7.82)  while  keeping  the  other  end  free.  Observe  that  the  attachment  line  on  sting 
is  all  in  the  upper  half  until  a bump  near  boattail-sting  juncture  is  developed  at  r=6.6 
to  drive  it  down.  At  t=7.8  a similar  bump  is  again  seen  near  the  juncture  but  to  pull 
the  string  up.  Note  that  when  such  a bump  is  formed  to  move  the  string  to  one  side,  the 
rear  part  of  the  string  is  still  gaining  its  magnitude  to  the  other  side.  Two  cycles  later 
the  amplitude  on  rear  sting  is  large  enough  to  interact  with  the  secondary  flow  region  as 
shown  on  the  top-right  corner  at  t=  10.6.  The  oscillation  of  skin-friction  lines  associated 
with  vortex  shedding,  including  the  primary/secondary  separation  and  attachment  lines 
are  shown  in  Figure  6.17  on  part  of  the  boattail  and  sting.  As  indicated  by  the  time 
history  of  crossflow  separation  in  Figure  5.27  and  the  ones  from  other  related  J cross- 
sections,  the  crossflow  secondary  flows  are  not  detected  for  x<  9.0  at  the  interested  times 
in  Figure  6.17.  Observe  the  merging  of  secondary  separation  and  attachment  lines  at 
7=9.5  on  top  half  of  the  plot  after  13.  With  the  primary  attachment  line  swinging 
to  the  lower  half  and  interacting  with  secondary  flow  for  r=9.9~  10.2,  the  secondary 
flow  region  on  the  top  half  shows  up  again.  The  process  reverses  itself  such  that  the 
upper  secondary  flow  fades  away  and  the  lower  one  re-appears.  Note  that  the  plot  at 
7=11.1  happens  to  be  very  close  to  the  mirror  image  of  the  one  at  7=10.2. 

The  crossflow  streamtraces  at  J = lll  (jc=10.85)  based  on  JMAX  = 131  grid  are 
presented  in  Figure  6.18  to  visualize  the  vortex  shedding  phenomena.  At  7=10.2  it  is 
clear  that  the  secondary  flow  does  not  exist  under  the  left,  smaller  vortex  which  stays 
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close  to  the  surface.  The  left  vortex  grows  in  size  and  pushes  the  right  vortex  away  from 
the  surface  (t=10.2~  10.6),  while  in  the  meantime,  a secondary  flow  region  is  formed 
under  the  left  vortex.  Although  the  size  of  left-hand-side  secondary  flow  is  practically 
invisible  at  this  moment,  it  can  be  identified  by  the  surface  skin-friction  lines  in  Figure 
6.17.  From  t=  10. 6 -10.9  it  appears  that  a right-hand-side  sub-vortex  center  is  left 
alone  near  the  surface  and  merges  into  the  seeondary  flow  to  form  a new  primary  vortex, 
while  the  previous  primary  vortex  is  dissipating  into  the  wake  region.  The  streamtraces 
shown  in  Figure  6.18  is  close  to  half  a cycle  by  observing  that  the  flow  pattern  at 
T=ll.l  is  similar  to  the  mirror  image  of  the  one  at  t=10.2.  Similarly  the  vortex 
shedding  patterns  at  J=115  (x=  11.88)  and  119  (x=  12.91)  are  presented  in  Figure  6.19 
and  6.20.  They  are  found  to  be  almost  the  same  as  the  ones  at  J=lll,  but  apparently 
with  a phase  lag.  Observe  that  the  flow  pattern  at  r=10.5  and  J = 115  is  very  close  to 
the  one  at  r=  10.3  and  J = 1 1 1 regarding  the  locations  of  primary  vortices,  indicating  that 
a lag  of  approximately  0.2  units  of  non-dimensional  time  is  present.  Similarly  the  flow 
pattern  at  J=119  clearly  falls  behind  the  one  at  J = 115  by  Ar=0. 2-0.3.  These  phase 
lags  suggest  that  the  vortex  shedding  on  boattail  is  not  parallel  to  the  x-axis,  but  rather 
inclined  with  a small  angle  to  the  sting.  The  vortical  flow  patterns  for  t=11.1  at  the 
cross-sections  J = 107  (x=9.83).  111,  115  and  119  are  shown  in  Figure  6.21.  Although 
the  shedding  is  not  seen  at  J=107,  the  streamtraces  at  the  other  3 cross-sections  clearly 
show  the  phase  lag  mentioned  above.  The  dense  streamtraces  within  the  left  primary 
vortex  at  J = 107  is  actually  an  unstable  limiting  eyele,  in  which  a stable  vortex  center  is 
not  shown.  Figure  6.22  presents  the  crossflow  streamtraces  at  t=  11.1  along  the  sting. 
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on  which  the  oblique  shedding  is  clearly  observed.  Toward  the  end  of  sting  the  shedding 
becomes  more  parallel.  As  a matter  of  fact  the  oblique  shedding  was  also  shown  in 
Ramberg’s  experiments  [19]  which  was  found  to  be  related  to  the  upstream  tip  vortices. 

Finally  in  Figure  6.23  the  oscillation  period  is  presented  along  the  sting  axis.  The 
period  is  obtained  from  the  time  history  of  crossflow  separation  points  based  on  the  latest 
simulated  oscillation  cycle.  For  those  locations  with  x<  10.5  the  periodic  oscillation  is 
also  observed  associated  with  mutual  interaction  between  two  primary  vortices,  which 
apparently  gives  a much  larger  period.  The  shedding  period  in  non-dimensional  time 
falls  between  1.75  to  1.9,  or  0.191  to  0.176  in  terms  of  the  Strouhal  number,  which  is 
very  comparable  to  the  value  0.182  given  by  Roshko’s  empirical  relation. 
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Figure  6.1  Evolution  of  crossflow  vortical  flow  at  J=54. 
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Figure  6.2  Development  of  crossflow  spiral  vortex  pattern  at  J=54. 
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Figure  6.3  Evolution  and  development  of  crossflow  vortical  flow  at  J=  19  (left) 

and  J=87  (right). 
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Figure  6.3  - continued. 
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Figure  6.4  Crossflow  streamtraces  at  J = 19  (top),  J=54  (middle)  and  J=87  (bottom) 

for  r=6.8  and  11.1. 
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Figure  6.5  Contour  plots  for  axial  velocity  component  (left)  and  crossflow  velocity 

magnitude  (right)  at  t=  11.1. 
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Figure  6.6  Contour  plots  for  duld^  (left)  and  pressure  coefficient  Cp  (right)  at  r=  11.1 
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Figure  6.7  Variation  of  cross-sectional  flow  topology  along  the  projectile 
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Figure  6.8  Time-sequence  of  surface  skin-friction  lines. 
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Figure  6.8-  continued 
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Figure  6.8  - continued. 
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Figure  6.9  Time-sequence  of  surface  Cp  distribution. 
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Figure  6.9  - continued. 
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Figure  6.9  - continued. 
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Figure  6. 10  Time-sequence  of  surface  temperature  (°R)  distribution. 
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Figure  6.10  - continued. 
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Figure  6. 10  - continued. 
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Figure  6. 1 1 Ogive  flow  development  - skin-friction  lines  (left) 

and  surface  Cp  distribution  (right). 
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Figure  6.11  - continued. 
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Figure  6. 12  Transition  of  longitudinal  flow  separation. 
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Figure  6.13  Convergence  of  surface  flow  patterns  - skin-friction  lines  (top),  Cp  (middle), 

and  temperature  (°R)  (bottom). 
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Figure  6.14  3-D  instantaneous  streamtraces  from  J=4  (x=0.086). 
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Figure  6.14  - continued. 
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Figure  6.14-  continued. 
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Figure  6. 15  3-D  streamtraces  from  J=35,  54,  68  and  77  at  t=  1 1 . 1 . 
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Figure  6.16  Oscillation  of  primary  attachment  line  near  0=180°. 
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Figure  6.16  - continued. 
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Figure  6.17  Oscillation  of  skin-friction  lines  on  boattail  (x<7.82)  and  sting  (x>7.82). 
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Figure  6.17  - continued. 
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Figure  6.17  - continued. 
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Figure  6.18  Alternating  vortex  shedding  at  J=I11  (x=10.85)  on  JMAX=131  grid. 
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Figure  6. 19  Alternating  vortex  shedding  at  J = 1 15  (x=  1 1 .88)  on  JMAX  =131  grid. 
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Figure  6.20  Alternating  vortex  shedding  at  J=  1 19  {x=  12.91)  on  JMAX=  131  grid. 
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Figure  6.21  Instantaneous  vortical  flow  pattern  on  sting  at  t=  1 1 . 1 . 


Figure  6.22  Variation  of  crossflow  vortex  shedding  pattern  on  sting  at  t=  1 1.1. 
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Figure  6.23  Oscillation  period  computed  on  projectile  sting. 


CHAPTER  VII 

SUMMARY  AND  CONCLUDING  REMARKS 


The  thin-layer  Navier-Stokes  simulation  of  an  impulsively-started  Mach  0.96 
turbulent  flow  past  a SCOBT  projectile  with  sting  at  60-degree  angle-of-attack  has  been 
conducted  by  implementing  the  original  Baldwin-Lomax  turbulence  model  and  a modified 
Baldwin-Lomax  model  proposed  by  Degani  and  Schiff.  The  original  Baldwin-Lomax 
model  is  an  algebraic  model  which  needs  not  to  solve  additional  equations  like  other  one- 
or  two-equation  models  do,  therefore  costs  little  extra  computer  time  which  is  especially 
advantageous  in  a large-scale  three-dimensional  flow  simulation.  It  is  known  that  many 
high  speed  flow  problems  can  be  accurately  simulated  by  the  thin-layer  Navier-Stokes 
computation  with  the  use  of  original  Baldwin-Lomax  turbulence  model,  provided  that  the 
flow  separation  is  not  very  strong.  In  the  high-a  vortical  flow  simulation  involving  large 
crossflow  separations,  however,  the  original  Baldwin-Lomax  model  was  found  to  over- 
predict the  eddy  viscosity  in  the  separated  region,  which  can  suppress  or  distort  the 
secondary  flow  phenomena  near  the  surface.  The  idea  behind  the  modified  Baldwin- 
Lomax  model  is  to  make  the  model  capable  of  differentiating  the  leeward  side  vortical 
flow  structure  by  properly  determining  a turbulent  mixing  length  within  the  turbulent 
boundary  layer,  therefore  effectively  reduce  the  eddy  viscosity  to  a reasonable  magnitude 
in  the  separated  region.  The  fundamental  formulations  in  the  original  Baldwin-Lomax 
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model  are  retained  such  that  the  vortical  flow  physics  and  phenomena  can  be  realistically 
modelled  without  requiring  further  computer  resources  or  tedious  coding  procedure.  The 
improvement  of  the  modified  Baldwin-Lomax  turbulence  model  has  been  demonstrated 
through  a number  of  subsonic  and  supersonic  flow  simulations  with  various  incidences 
in  the  past  decade.  The  similar  investigations  for  transonic  high-a  flows,  however,  were 
very  limited. 

Although  there  are  no  experimental  measurements  available  to  justify  the 
numerical  accuracy  of  the  current  computational  work,  a study  to  the  effects  of  original 
Baldwin-Lomax  turbulence  model  and  the  modified  model  upon  the  general  vortical  flow 
characteristics  has  revealed  that  the  modified  model  is  indeed  capable  of  bringing  up  the 
complex  secondary  flow  phenomena.  The  simulated  flowfield  by  modified  Baldwin- 
Lomax  model  has  been  shown  to  be  quite  comparable  in  quality  to  the  results  from  other 
high-a  projectile  flow  studies.  Some  of  the  important  findings  are  summarized  as 
follows. 

(1)  The  amount  of  computed  eddy  viscosity  is  significantly  reduced  in  the 
separated  region  by  using  the  modified  Baldwin-Lomax  turbulence  model.  The 
secondary  crossflow  separations  are  simulated  over  the  projectile  surface  which  are 
mostly  absent  from  the  original  Baldwin-Lomax  model  computation.  The  diminishing 
of  strong  diffusive  process  in  the  separated  region  has  altered  the  general  flow 
characteristics  in  the  boundary  layer.  The  flow  in  the  primary  vortex  region  away  from 
the  surface,  on  the  other  hand,  is  almost  unaffected  despite  of  the  substantial  reduction 
of  eddy  viscosity  magnitude  in  the  modified  model  computation. 
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(2)  Although  the  modified  Baldwin-Lomax  model  has  predicted  a more  complex 
surface  flow  pattern,  the  resulting  surface  pressure  distribution  is  found  to  be  very 
comparable  to  the  one  obtained  from  original  model  computation  for  an  established  flow. 
The  simulated  secondary  flows  are  apparently  too  weak  to  induce  significant  surface 
pressure  difference.  The  viscous  force  does  change  with  a different  surface  skin-friction 
lines  distribution  for  modified  Baldwin-Lomax  model  results,  but  consists  of  barely  1 % 
of  the  pressure  force  in  magnitude.  Accordingly,  the  overall  lift  and  drag  forces  are  less 
affected.  A fluctuating  side  force  time  history  is  observed  due  to  the  development  of 
alternating  vortex  shedding  on  sting.  The  oscillation  is  about  a zero  mean,  with  an 
amplitude  of  4%  to  5%  of  the  lift  forces. 

(The  followings  are  based  on  modified  Baldwin-Lomax  model  only.) 

(3)  The  primary  crossflow  separation  starts  to  develop  from  the  leeward 
symmetric  plane  and  moves  upstream  along  the  circumferential  direction.  The  onset  time 
for  primary  separation  is  found  to  be  linearly  dependent  to  the  local  projectile  diameter. 
The  crossflow  secondary  separation  is  brought  into  existence  when  enough  adverse 
pressure  gradient  is  developed  under  the  primary  vortex  structure.  The  onset  time  for 
secondary  separation  is  also  linearly  dependent  to  the  diameter,  with  a slight  delay  for 
those  cross-sections  on  boattail  compared  to  the  ones  on  ogive  having  the  same 
diameters. 

(4)  The  evolution  of  crossflow  topology  on  the  cylinder,  boattail  and  sting  sections 
involves  transition  of  a stable  spiral  vortex  to  an  unstable  one  through  the  appearance  of 
a stable  limiting  cycle.  Only  one  saddle  point  is  located  on  the  leeward  symmetric  plane 
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to  separate  the  upward  and  downward  flows.  The  evolution  on  the  ogive  section,  on  the 
other  hand,  consists  of  two  saddle  points  on  the  leeward  symmetric  plane  and  one 
unstable  spiral  vortex.  The  two-saddle-point  structure  is  found  to  be  related  to  the 
surface  curvature  in  conjunction  with  the  local  flow  direction. 

(5)  The  circumferential  velocity  component  dominates  the  surface  flow  pattern  on 
the  windward  side  before  crossflow  primary  separation.  Accordingly,  the  surface  skin- 
friction  lines  are  almost  perpendicular  to  the  primary  separation  line.  The  induced 
crossflow  secondary  flows  in  the  early  flow  development  are  only  found  on  smaller 
cross-sections  on  ogive  and  boattail-sting,  which  gradually  propagate  onto  the  cylindrical 
section  and  merge  together  as  time  grows.  A high  speed  Jet  stream  having  large  axial 
velocity  (almost  twice  of  the  free  stream  axial  velocity)  is  induced  above  the  leeward 
symmetric  plane  over  ogive  and  cylinder,  which  provides  an  effective  "channel"  to  bring 
massive  longitudinal  flow  from  the  ogive  to  the  boattail  and  sting. 

(6)  The  particular  boattail  geometry  has  caused  a longitudinally  reverse  flow  right 
after  the  cylinder-boattail  Juncture  in  the  crossflow  separated  region.  Behind  that  the 
secondary  flow  structure  collapses  due  to  a relatively  weak  circumferential  flow  action 
until  X » 10  on  sting  where  the  circumferential  component  becomes  strong  enough  to 
induce  the  adverse  pressure  gradient  and  the  secondary  flow  separation  re-appears. 

(7)  The  vortex  shedding  is  simulated  only  on  the  last  2/3  of  sting  section,  which 
is  believed  to  be  a result  of  absolute  instabilities.  The  Strouhal  number  is  in  excellent 
agreement  to  a classical  empirical  value  of  0.182.  The  shedding  is  found  to  be  slightly 
oblique  near  the  boattail  and  gradually  becomes  parallel  near  the  end  of  sting. 


167 


The  investigation  in  this  study  has  revealed  some  important  information  to  the 
high-a  projectile  aerodynamics.  The  modified  Baldwin-Lomax  turbulence  model  is 
justified  to  be  capable  of  representing  the  high  speed  turbulent  vortical  flow  physics  and 
phenomena.  The  choice  of  cut-off  constant  needs  to  be  properly  resolved  in  order  to 
have  more  quantitatively  accurate  solution.  In  order  to  gain  a more  comprehensive 
understanding  to  the  complex  three-dimensional  vortical  flow  phenomena,  it  is  desired 
to  have  further  computational  investigations  to  the  flow  cases  with  various  angles  of 
attack,  Mach  numbers  and  Reynolds  numbers. 


APPENDIX 

JACOBIAN  MATRICES 


The  Jacobian  matrices  A,  B and  C have  the  following  form 
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Here  k = and  f to  represent  A,  B and  C,  respectively.  Note  that  A:,  = 0 for  a fixed 
grid  case. 
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where 


The  Jacobian  matrix  M has  the  following  form 
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